R Markdown: The Definitive Guide https://bookdown.org/yihui/rmarkdown/
Tools for Reproducible Research http://kbroman.org/Tools4RR/
options(scipen = 999) ## disable scientific notation 4.835312e-04 --> 0.0004835312
library(tidyverse)
#> ✔ ggplot2 2.2.1.9000 ✔ purrr 0.2.4
#> ✔ tibble 1.4.2 ✔ dplyr 0.7.4
#> ✔ tidyr 0.8.0 ✔ stringr 1.3.0
#> ✔ readr 1.1.1 ✔ forcats 0.3.0
library(readxl) #read excel
library(stargazer) #nice and informative tables
library(compareGroups) # nice table one
library(summarytools) # summary table for exploratory analysis
library(labelVector) # labeling variables
library(PoEdata)#for PoE datasets
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot. You can use include=FALSE to exclude everything in a chunk. If you only want to suppress messages, use message=FALSE instead. We can mask the result by using results = FALSE. Use command-shift-m to insert pipe operator %>%.
#rename several variables
#dplyr::rename(.data, NEW.NAME1 = old.name1, NEW.NAME2 = old.name2)
## An 'if' statement takes the form:
## if(logical boolean){function to perform if the boolean is TRUE}
x <- 100000000
if(x > 10000){
"x is a really big number"
}
## [1] "x is a really big number"
if(x < 10000){
"x isn't that big"
} ## Nothing happens here because x < 10000 evaluates to false
## We can use else statements to deal with instances where the boolean
## evaluates to FALSE
## Note that using multiple if/else statements requires that you wrap the whole thing
## in curly brackets, {}:
x <- 1
if(x > 10000){
"x is a really big number"
} else{
"x isn't that big"
}
## [1] "x isn't that big"
#######################################################
## Note about the modulus function: what it does is divides the number
## in front of it by the number after it, and gives you the remainder.
## It's really useful to use inside of loops because you can use it to
## print the iteration you're on.
## In other words, if you have a loop that has 100,000 iterations,
## maybe you want to keep tabs on how long it's taking and where it is,
## but don't want to print every single number from 1 to 100,000.
## Inserting this code into the loop will print only every 100th iteration:
for(i in 1:500){
if((i %% 100) == 0){
print(i)
}
}
## [1] 100
## [1] 200
## [1] 300
## [1] 400
## [1] 500
####################################################
## For a loop that does 10,000 iterations:
start <- proc.time()[3] ## this captures and stores the time before the loop starts
holder <- c()
for(i in 1:10000){
holder[i] <- rnorm(1, mean = 0, sd = 1)
# here we are drawing one observation from a normal distribution
# with mean zero and standard deviation 1
}
end <- proc.time()[3]
## How long did that take?
loop10000 <- end - start
loop10000
## elapsed
## 0.035
## What if we had vectorized our code and drawn all 10,000 random numbers first?
start <- proc.time()[3]
draws <- rnorm(10000, mean = 0, sd = 1)
end <- proc.time()[3]
vect10000 <- end - start
loop10000
## elapsed
## 0.035
vect10000
## elapsed
## 0.001
## The vectorized code was way, way faster!
####################################################
## There's one other way, and it's a really cool R-hack.
## When you write a function, you can store objects using
## <<- instead of just <-
## Using the double arrow stores the object to your workspace,
## so that it's accessible after the function is complete
debugging.example <- function(x){
holder <- rbeta(1000, 1, 10)
holder <- sort(holder, decreasing = T)
holder <- matrix(holder[1:50], 25, 2)
# put the <<- here to store holder
holder <<- matrix(holder[1:50], 25, 2)
t(holder) %*% x
}
#debugging.example(test.matrix) ## Now it works
head(holder)
## [1] 0.38303320 -0.49538108 -0.14127816 0.01735916 0.64212009 -0.96262673
# Another trick to debug functions is to store objects in it with
# '<<-' instead of '<-'. If you only use '<-' then the only thing
# left in your workspace after the function runs will be the output
# from the final line of the function. Using '<<-' will force that
# object into your workspace, even if it's not the last line in
# the function
####################################################
# Note that with lists you can use lapply() or sapply()
# lapply() will output your result as a list
# sapply() will (usually) output your result as a vector
####################################################
## Debugging functions
## The browser() function is super super helpful if you have a really
## complicated function that has a bug in the code.
## Now let's see how we could use the browser() function to debug.
## Note: you can also put browser() functions into for/while loops
debugging.example <- function(x){
holder <- rbeta(1000, 1, 10)
browser()
holder <- sort(holder, decreasing = T)
browser()
holder <- matrix(holder[1:50], 25, 2)
browser()
holder %*% x
}
test.matrix <- matrix(rnorm(100), 25, 4)
debugging.example(test.matrix)
## this tells us that the error doesn't occur until after the 3rd browser() statement
library(parallel)
library(doParallel)
## Loading required package: foreach
##
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
##
## accumulate, when
## Loading required package: iterators
library(foreach)
#Starting a cluster
numCores <- detectCores()
numCores
## [1] 8
cl <- makeCluster(numCores-1)
#Close the cluster
stopCluster(cl)
#Introduction to summarytools
#https://cran.r-project.org/web/packages/summarytools/vignettes/Introduction.html
#install.packages('devtools')
#library(devtools)
#install_github('dcomtois/summarytools')
#view(dfSummary(iris)) #open a new browser to view the summary
#graph.magnif-> graph size
print(dfSummary(tobacco, style = 'grid', plain.ascii = FALSE, graph.magnif = 4),
method = 'render', headings = FALSE)
| No | Variable | Stats / Values | Freqs (% of Valid) | Graph | Valid | Missing |
|---|---|---|---|---|---|---|
| 1 | gender [factor] | 1. F 2. M | 489 (50.0%) 489 (50.0%) | 978 (97.8%) | 22 (2.2%) | |
| 2 | age [numeric] | mean (sd) : 49.6 (18.29) min < med < max : 18 < 50 < 80 IQR (CV) : 32 (0.37) | 63 distinct values | 975 (97.5%) | 25 (2.5%) | |
| 3 | age.gr [factor] | 1. 18-34 2. 35-50 3. 51-70 4. 71 + | 258 (26.5%) 241 (24.7%) 317 (32.5%) 159 (16.3%) | 975 (97.5%) | 25 (2.5%) | |
| 4 | BMI [numeric] | mean (sd) : 25.73 (4.49) min < med < max : 8.83 < 25.62 < 39.44 IQR (CV) : 5.72 (0.17) | 974 distinct values | 974 (97.4%) | 26 (2.6%) | |
| 5 | smoker [factor] | 1. Yes 2. No | 298 (29.8%) 702 (70.2%) | 1000 (100%) | 0 (0%) | |
| 6 | cigs.per.day [numeric] | mean (sd) : 6.78 (11.88) min < med < max : 0 < 0 < 40 IQR (CV) : 11 (1.75) | 37 distinct values | 965 (96.5%) | 35 (3.5%) | |
| 7 | diseased [factor] | 1. Yes 2. No | 224 (22.4%) 776 (77.6%) | 1000 (100%) | 0 (0%) | |
| 8 | disease [character] | 1. Hypertension 2. Cancer 3. Cholesterol 4. Heart 5. Pulmonary 6. Musculoskeletal 7. Diabetes 8. Hearing 9. Digestive 10. Hypotension [ 3 others ] | 36 (16.2%) 34 (15.3%) 21 (9.5%) 20 (9.0%) 20 (9.0%) 19 (8.6%) 14 (6.3%) 14 (6.3%) 12 (5.4%) 11 (5.0%) 21 (9.5%) | 222 (22.2%) | 778 (77.8%) | |
| 9 | samp.wgts [numeric] | mean (sd) : 1 (0.08) min < med < max : 0.86 < 1.04 < 1.06 IQR (CV) : 0.19 (0.08) | 0.86! : 267 (26.7%) 1.04! : 249 (24.9%) 1.05! : 324 (32.4%) 1.06! : 160 (16.0%) ! rounded | 1000 (100%) | 0 (0%) |
Generated by summarytools 0.8.9 (R version 3.5.2)
2019-01-07
#http://lutgw1.lunet.edu/usr/lib64/R/library/Hmisc/tests/examples.Rmd
#http://hbiostat.org/R/Hmisc/examples.html
library(Hmisc)
d<-describe(tobacco)
html(d, size=90, scroll=F)
| n | missing | distinct |
|---|---|---|
| 978 | 22 | 2 |
Value F M Frequency 489 489 Proportion 0.5 0.5
| n | missing | distinct | Info | Mean | Gmd | .05 | .10 | .25 | .50 | .75 | .90 | .95 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 975 | 25 | 63 | 1 | 49.6 | 21.12 | 21.0 | 24.0 | 34.0 | 50.0 | 66.0 | 74.0 | 77.3 |
| n | missing | distinct |
|---|---|---|
| 975 | 25 | 4 |
Value 18-34 35-50 51-70 71 + Frequency 258 241 317 159 Proportion 0.265 0.247 0.325 0.163
| n | missing | distinct | Info | Mean | Gmd | .05 | .10 | .25 | .50 | .75 | .90 | .95 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 974 | 26 | 974 | 1 | 25.73 | 5.039 | 18.47 | 20.24 | 22.93 | 25.62 | 28.65 | 31.55 | 33.11 |
| lowest : | 8.825777 | 9.008713 | 10.346131 | 12.984157 | 14.179153 |
| highest: | 38.122074 | 38.372244 | 38.637694 | 39.207753 | 39.439012 |
| n | missing | distinct |
|---|---|---|
| 1000 | 0 | 2 |
Value Yes No Frequency 298 702 Proportion 0.298 0.702
| n | missing | distinct | Info | Mean | Gmd | .05 | .10 | .25 | .50 | .75 | .90 | .95 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 965 | 35 | 37 | 0.659 | 6.782 | 10.6 | 0 | 0 | 0 | 0 | 11 | 29 | 35 |
| n | missing | distinct |
|---|---|---|
| 1000 | 0 | 2 |
Value Yes No Frequency 224 776 Proportion 0.224 0.776
| n | missing | distinct |
|---|---|---|
| 222 | 778 | 13 |
| n | missing | distinct | Info | Mean | Gmd |
|---|---|---|---|---|---|
| 1000 | 0 | 4 | 0.927 | 1 | 0.07774 |
Value 0.8614232 1.0441767 1.0493827 1.0625000 Frequency 267 249 324 160 Proportion 0.267 0.249 0.324 0.160
#compareGroups 4.0: Descriptives by groups
#https://cran.r-project.org/web/packages/compareGroups/vignettes/compareGroups_vignette.html
data(predimed)
predimed$tmain <- with(predimed, Surv(toevent, event == "Yes"))
Hmisc::label(predimed$tmain) <- "AMI, stroke, or CV Death"
# https://cran.r-project.org/web/packages/labelVector/vignettes/labelVector.html
labelVector::get_label(predimed)
## [1] "Intervention group" "Sex"
## [3] "Age" "Smoking"
## [5] "Body mass index" "Waist circumference"
## [7] "Waist-to-height ratio" "Hypertension"
## [9] "Type-2 diabetes" "Dyslipidemia"
## [11] "Family history of premature CHD" "Hormone-replacement therapy"
## [13] "MeDiet Adherence score" "follow-up to main event (years)"
## [15] "AMI, stroke, or CV Death" "AMI, stroke, or CV Death"
#set label
x <- 1:10
x <- set_label(x, "some integers")
x
## some integers
## [1] 1 2 3 4 5 6 7 8 9 10
#To remove the variable toevent and event from the analysis:
compareGroups(group ~ . - toevent - event, data = predimed)
##
##
## -------- Summary of results by groups of 'Intervention group'---------
##
##
## var N p.value method
## 1 Sex 6324 <0.001** categorical
## 2 Age 6324 0.003** continuous normal
## 3 Smoking 6324 0.444 categorical
## 4 Body mass index 6324 <0.001** continuous normal
## 5 Waist circumference 6324 0.045** continuous normal
## 6 Waist-to-height ratio 6324 <0.001** continuous normal
## 7 Hypertension 6324 0.249 categorical
## 8 Type-2 diabetes 6324 0.017** categorical
## 9 Dyslipidemia 6324 0.423 categorical
## 10 Family history of premature CHD 6324 0.581 categorical
## 11 Hormone-replacement therapy 5661 0.850 categorical
## 12 MeDiet Adherence score 6324 <0.001** continuous normal
## 13 AMI, stroke, or CV Death 6324 0.011** Surv [Tmax=4.79]
## selection
## 1 ALL
## 2 ALL
## 3 ALL
## 4 ALL
## 5 ALL
## 6 ALL
## 7 ALL
## 8 ALL
## 9 ALL
## 10 ALL
## 11 ALL
## 12 ALL
## 13 ALL
## -----
## Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1
No filters have been used (e.g., selecting only treated patients); therefore, the selection column lists ALL (for all variables).
#To perform the analysis in a subset of participants (e.g., “female” participants):
compareGroups(group ~ age + smoke + waist + hormo, data = predimed,
subset = sex == "Female")
##
##
## -------- Summary of results by groups of 'group'---------
##
##
## var N p.value method
## 1 Age 3645 0.056* continuous normal
## 2 Smoking 3645 0.907 categorical
## 3 Waist circumference 3645 0.016** continuous normal
## 4 Hormone-replacement therapy 3459 0.898 categorical
## selection
## 1 sex == "Female"
## 2 sex == "Female"
## 3 sex == "Female"
## 4 sex == "Female"
## -----
## Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1
#To subset specific variable/s (e.g., hormo and waist):
compareGroups(group ~ age + sex + smoke + waist + hormo, data = predimed,
selec = list(hormo = sex == "Female", waist = waist > 20))
##
##
## -------- Summary of results by groups of 'Intervention group'---------
##
##
## var N p.value method
## 1 Age 6324 0.003** continuous normal
## 2 Sex 6324 <0.001** categorical
## 3 Smoking 6324 0.444 categorical
## 4 Waist circumference 6324 0.045** continuous normal
## 5 Hormone-replacement therapy 3459 0.898 categorical
## selection
## 1 ALL
## 2 ALL
## 3 ALL
## 4 waist > 20
## 5 sex == "Female"
## -----
## Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1
#Combinations are also allowed, e.g.:
compareGroups(group ~ age + smoke + waist + hormo, data = predimed,
selec = list(waist = !is.na(hormo)), subset = sex == "Female")
##
##
## -------- Summary of results by groups of 'group'---------
##
##
## var N p.value method
## 1 Age 3645 0.056* continuous normal
## 2 Smoking 3645 0.907 categorical
## 3 Waist circumference 3459 0.007** continuous normal
## 4 Hormone-replacement therapy 3459 0.898 categorical
## selection
## 1 sex == "Female"
## 2 sex == "Female"
## 3 (sex == "Female") & (!is.na(hormo))
## 4 sex == "Female"
## -----
## Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1
#A variable can appear twice in the formula, e.g.: (NOT SHOW)
#To change default options, e.g., “waist” used as non-normal distributed:
compareGroups(group ~ age + smoke + waist + hormo, data = predimed,
method = c(waist = 2))
## Warning in cor.test.default(x, as.integer(y), method = "spearman"): Cannot
## compute exact p-value with ties
##
##
## -------- Summary of results by groups of 'Intervention group'---------
##
##
## var N p.value method selection
## 1 Age 6324 0.003** continuous normal ALL
## 2 Smoking 6324 0.444 categorical ALL
## 3 Waist circumference 6324 0.085* continuous non-normal ALL
## 4 Hormone-replacement therapy 5661 0.850 categorical ALL
## -----
## Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1
#If the method argument is stated as NA for a variable, then a Shapiro-Wilk test for normality is used to decide if the variable is normal or non-normal distributed. To change the significance threshold:
compareGroups(group ~ age + smoke + waist + hormo, data = predimed,
method = c(waist = NA), alpha = 0.01)
## Warning in cor.test.default(x, as.integer(y), method = "spearman"): Cannot
## compute exact p-value with ties
##
##
## -------- Summary of results by groups of 'Intervention group'---------
##
##
## var N p.value method selection
## 1 Age 6324 0.003** continuous normal ALL
## 2 Smoking 6324 0.444 categorical ALL
## 3 Waist circumference 6324 0.085* continuous non-normal ALL
## 4 Hormone-replacement therapy 5661 0.850 categorical ALL
## -----
## Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1
#Shapiro-Wilk test, stating the cutpoint at 0.01 level
#1: forces analysis as normal-distributed
#2: forces analysis as continuous non-normal
#3: forces analysis as categorical
#NA: performs a Shapiro-Wilks test to decide between normal or non-normal
#All non factor variables are considered as continuous. Exception is made (by default) for those that have fewer than 5 different values.
#The object from compareGroups can later be updated. For example:
res <- compareGroups(group ~ age + sex + smoke + waist + hormo,
data = predimed)
res <- update(res, . ~ . - sex + bmi + toevent, subset = sex ==
"Female", method = c(waist = 2, tovent = 2), selec = list(bmi = !is.na(hormo)))
## Warning in compareGroups.fit(X = X, y = y, include.label = include.label, :
## variables tovent specified in 'method' not found
## Warning in compareGroups.fit(X = X, y = y, include.label = include.label, :
## Cannot compute exact p-value with ties
res
##
##
## -------- Summary of results by groups of 'group'---------
##
##
## var N p.value method
## 1 Age 3645 0.056* continuous normal
## 2 Smoking 3645 0.907 categorical
## 3 Waist circumference 3645 0.037** continuous non-normal
## 4 Hormone-replacement therapy 3459 0.898 categorical
## 5 Body mass index 3459 0.002** continuous normal
## 6 follow-up to main event (years) 3645 <0.001** continuous normal
## selection
## 1 sex == "Female"
## 2 sex == "Female"
## 3 sex == "Female"
## 4 sex == "Female"
## 5 (sex == "Female") & (!is.na(hormo))
## 6 sex == "Female"
## -----
## Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1
#the Odds Ratio (OR) can be printed in the final table. If the response variable is time-to-event (see Section 3.1), the Hazard Ratio (HR) can be printed instead.
res2 <- compareGroups(htn ~ age + sex + bmi + smoke, data = predimed,
ref = c(smoke = 1, sex = 2))
createTable(res2, show.ratio = TRUE)
##
## --------Summary descriptives table by 'Hypertension'---------
##
## ___________________________________________________________________________
## No Yes OR p.ratio p.overall
## N=1089 N=5235
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## Age 65.9 (6.19) 67.2 (6.15) 1.04 [1.03;1.05] <0.001 <0.001
## Sex: <0.001
## Male 595 (54.6%) 2084 (39.8%) 0.55 [0.48;0.63] 0.000
## Female 494 (45.4%) 3151 (60.2%) Ref. Ref.
## Body mass index 28.9 (3.69) 30.2 (3.80) 1.10 [1.08;1.12] <0.001 <0.001
## Smoking: <0.001
## Never 536 (49.2%) 3356 (64.1%) Ref. Ref.
## Current 233 (21.4%) 625 (11.9%) 0.43 [0.36;0.51] 0.000
## Former 320 (29.4%) 1254 (24.0%) 0.63 [0.54;0.73] <0.001
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
#The createTable function, applied to an object of compareGroups class, returns tables with descriptives that can be displayed on-screen or exported to CSV, LaTeX, HTML, Word or Excel.
res <- compareGroups(group ~ age + sex + smoke + waist + hormo,
data = predimed, selec = list(hormo = sex == "Female"))
restab <- createTable(res)
#the option “descr” prints descriptive tables.
#digits: The number of digits that appear in the results can be changed
print(restab, which.table = "descr")
##
## --------Summary descriptives table by 'Intervention group'---------
##
## ________________________________________________________________________________
## Control MedDiet + Nuts MedDiet + VOO p.overall
## N=2042 N=2100 N=2182
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## Age 67.3 (6.28) 66.7 (6.02) 67.0 (6.21) 0.003
## Sex: <0.001
## Male 812 (39.8%) 968 (46.1%) 899 (41.2%)
## Female 1230 (60.2%) 1132 (53.9%) 1283 (58.8%)
## Smoking: 0.444
## Never 1282 (62.8%) 1259 (60.0%) 1351 (61.9%)
## Current 270 (13.2%) 296 (14.1%) 292 (13.4%)
## Former 490 (24.0%) 545 (26.0%) 539 (24.7%)
## Waist circumference 101 (10.8) 100 (10.6) 100 (10.4) 0.045
## Hormone-replacement therapy: 0.898
## No 1143 (97.4%) 1036 (97.2%) 1183 (97.0%)
## Yes 31 (2.64%) 30 (2.81%) 36 (2.95%)
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
#show.all: If show.all argument is set to TRUE a column is displayed with descriptives for all data:
createTable(res, show.all = TRUE)
##
## --------Summary descriptives table by 'Intervention group'---------
##
## _____________________________________________________________________________________________
## [ALL] Control MedDiet + Nuts MedDiet + VOO p.overall
## N=6324 N=2042 N=2100 N=2182
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## Age 67.0 (6.17) 67.3 (6.28) 66.7 (6.02) 67.0 (6.21) 0.003
## Sex: <0.001
## Male 2679 (42.4%) 812 (39.8%) 968 (46.1%) 899 (41.2%)
## Female 3645 (57.6%) 1230 (60.2%) 1132 (53.9%) 1283 (58.8%)
## Smoking: 0.444
## Never 3892 (61.5%) 1282 (62.8%) 1259 (60.0%) 1351 (61.9%)
## Current 858 (13.6%) 270 (13.2%) 296 (14.1%) 292 (13.4%)
## Former 1574 (24.9%) 490 (24.0%) 545 (26.0%) 539 (24.7%)
## Waist circumference 100 (10.6) 101 (10.8) 100 (10.6) 100 (10.4) 0.045
## Hormone-replacement therapy: 0.898
## No 3362 (97.2%) 1143 (97.4%) 1036 (97.2%) 1183 (97.0%)
## Yes 97 (2.80%) 31 (2.64%) 30 (2.81%) 36 (2.95%)
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
#Tables made with the same response variable can be combined by row:
restab1 <- createTable(compareGroups(group ~ age + sex, data = predimed))
restab2 <- createTable(compareGroups(group ~ bmi + smoke, data = predimed))
rbind(`Non-modifiable risk factors` = restab1, `Modifiable risk factors` = restab2)
##
## --------Summary descriptives table by 'Intervention group'---------
##
## _______________________________________________________________________
## Control MedDiet + Nuts MedDiet + VOO p.overall
## N=2042 N=2100 N=2182
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## Non-modifiable risk factors:
## Age 67.3 (6.28) 66.7 (6.02) 67.0 (6.21) 0.003
## Sex: <0.001
## Male 812 (39.8%) 968 (46.1%) 899 (41.2%)
## Female 1230 (60.2%) 1132 (53.9%) 1283 (58.8%)
## Modifiable risk factors:
## Body mass index 30.3 (3.96) 29.7 (3.77) 29.9 (3.71) <0.001
## Smoking: 0.444
## Never 1282 (62.8%) 1259 (60.0%) 1351 (61.9%)
## Current 270 (13.2%) 296 (14.1%) 292 (13.4%)
## Former 490 (24.0%) 545 (26.0%) 539 (24.7%)
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
#Columns from tables built with the same explanatory and response variables but done with a different subset (i.e. “Male” and “Female”, strata) can be combined:
# first build the table with descriptives by groups:
# include.miss = T shows missing value category
res <- compareGroups(group ~ . - sex, include.miss = T, predimed)
restab <- createTable(res, hide.no = "no")
# and then apply the strataTable function on the table:
strataTable(restab, "sex")
##
## --------Summary descriptives table ---------
##
## _______________________________________________________________________________________________________________________________________
## Male Female
## __________________________________________________ ___________________________________________________
## Control MedDiet + Nuts MedDiet + VOO p.overall Control MedDiet + Nuts MedDiet + VOO p.overall
## N=812 N=968 N=899 N=1230 N=1132 N=1283
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## Age 66.4 (6.62) 65.8 (6.40) 66.1 (6.61) 0.215 68.0 (5.96) 67.4 (5.57) 67.7 (5.84) 0.056
## Smoking: 0.851 0.907
## Never 205 (25.2%) 266 (27.5%) 236 (26.3%) 1077 (87.6%) 993 (87.7%) 1115 (86.9%)
## Current 204 (25.1%) 242 (25.0%) 221 (24.6%) 66 (5.37%) 54 (4.77%) 71 (5.53%)
## Former 403 (49.6%) 460 (47.5%) 442 (49.2%) 87 (7.07%) 85 (7.51%) 97 (7.56%)
## Body mass index 29.6 (3.45) 29.1 (3.28) 29.2 (3.28) 0.018 30.8 (4.20) 30.2 (4.08) 30.4 (3.91) 0.002
## Waist circumference 104 (9.82) 103 (9.36) 103 (9.65) 0.289 99.0 (11.0) 97.8 (11.0) 98.0 (10.5) 0.016
## Waist-to-height ratio 0.62 (0.06) 0.62 (0.06) 0.62 (0.06) 0.191 0.64 (0.07) 0.63 (0.07) 0.63 (0.07) 0.002
## Hypertension 649 (79.9%) 753 (77.8%) 682 (75.9%) 0.130 1062 (86.3%) 985 (87.0%) 1104 (86.0%) 0.780
## Type-2 diabetes 430 (53.0%) 496 (51.2%) 486 (54.1%) 0.468 540 (43.9%) 454 (40.1%) 596 (46.5%) 0.007
## Dyslipidemia 531 (65.4%) 653 (67.5%) 589 (65.5%) 0.575 948 (77.1%) 886 (78.3%) 971 (75.7%) 0.319
## Family history of premature CHD 135 (16.6%) 171 (17.7%) 156 (17.4%) 0.841 327 (26.6%) 289 (25.5%) 351 (27.4%) 0.596
## Hormone-replacement therapy: 0.903 0.689
## No 668 (82.3%) 799 (82.5%) 735 (81.8%) 1143 (92.9%) 1036 (91.5%) 1183 (92.2%)
## Yes 0 (0.00%) 0 (0.00%) 0 (0.00%) 31 (2.52%) 30 (2.65%) 36 (2.81%)
## 'Missing' 144 (17.7%) 169 (17.5%) 164 (18.2%) 56 (4.55%) 66 (5.83%) 64 (4.99%)
## MeDiet Adherence score 8.57 (1.94) 8.86 (1.96) 8.91 (2.02) 0.001 8.36 (1.94) 8.77 (1.86) 8.68 (1.92) <0.001
## follow-up to main event (years) 4.05 (1.78) 4.38 (1.74) 4.53 (1.64) <0.001 4.12 (1.71) 4.26 (1.67) 4.72 (1.56) <0.001
## AMI, stroke, or CV Death 58 (7.14%) 41 (4.24%) 52 (5.78%) 0.029 39 (3.17%) 29 (2.56%) 33 (2.57%) 0.576
## AMI, stroke, or CV Death 8.48% 4.42% 5.71% 0.008 3.84% 2.66% 2.16% 0.269
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
export2md(strataTable(restab, "sex"), size = 8)
| Control | MedDiet + Nuts | MedDiet + VOO | p.overall | Control | MedDiet + Nuts | MedDiet + VOO | p.overall | |
|---|---|---|---|---|---|---|---|---|
| N=812 | N=968 | N=899 | N=1230 | N=1132 | N=1283 | |||
| Age | 66.4 (6.62) | 65.8 (6.40) | 66.1 (6.61) | 0.215 | 68.0 (5.96) | 67.4 (5.57) | 67.7 (5.84) | 0.056 |
| Smoking: | 0.851 | 0.907 | ||||||
| Never | 205 (25.2%) | 266 (27.5%) | 236 (26.3%) | 1077 (87.6%) | 993 (87.7%) | 1115 (86.9%) | ||
| Current | 204 (25.1%) | 242 (25.0%) | 221 (24.6%) | 66 (5.37%) | 54 (4.77%) | 71 (5.53%) | ||
| Former | 403 (49.6%) | 460 (47.5%) | 442 (49.2%) | 87 (7.07%) | 85 (7.51%) | 97 (7.56%) | ||
| Body mass index | 29.6 (3.45) | 29.1 (3.28) | 29.2 (3.28) | 0.018 | 30.8 (4.20) | 30.2 (4.08) | 30.4 (3.91) | 0.002 |
| Waist circumference | 104 (9.82) | 103 (9.36) | 103 (9.65) | 0.289 | 99.0 (11.0) | 97.8 (11.0) | 98.0 (10.5) | 0.016 |
| Waist-to-height ratio | 0.62 (0.06) | 0.62 (0.06) | 0.62 (0.06) | 0.191 | 0.64 (0.07) | 0.63 (0.07) | 0.63 (0.07) | 0.002 |
| Hypertension | 649 (79.9%) | 753 (77.8%) | 682 (75.9%) | 0.130 | 1062 (86.3%) | 985 (87.0%) | 1104 (86.0%) | 0.780 |
| Type-2 diabetes | 430 (53.0%) | 496 (51.2%) | 486 (54.1%) | 0.468 | 540 (43.9%) | 454 (40.1%) | 596 (46.5%) | 0.007 |
| Dyslipidemia | 531 (65.4%) | 653 (67.5%) | 589 (65.5%) | 0.575 | 948 (77.1%) | 886 (78.3%) | 971 (75.7%) | 0.319 |
| Family history of premature CHD | 135 (16.6%) | 171 (17.7%) | 156 (17.4%) | 0.841 | 327 (26.6%) | 289 (25.5%) | 351 (27.4%) | 0.596 |
| Hormone-replacement therapy: | 0.903 | 0.689 | ||||||
| No | 668 (82.3%) | 799 (82.5%) | 735 (81.8%) | 1143 (92.9%) | 1036 (91.5%) | 1183 (92.2%) | ||
| Yes | 0 (0.00%) | 0 (0.00%) | 0 (0.00%) | 31 (2.52%) | 30 (2.65%) | 36 (2.81%) | ||
| ‘Missing’ | 144 (17.7%) | 169 (17.5%) | 164 (18.2%) | 56 (4.55%) | 66 (5.83%) | 64 (4.99%) | ||
| MeDiet Adherence score | 8.57 (1.94) | 8.86 (1.96) | 8.91 (2.02) | 0.001 | 8.36 (1.94) | 8.77 (1.86) | 8.68 (1.92) | <0.001 |
| follow-up to main event (years) | 4.05 (1.78) | 4.38 (1.74) | 4.53 (1.64) | <0.001 | 4.12 (1.71) | 4.26 (1.67) | 4.72 (1.56) | <0.001 |
| AMI, stroke, or CV Death | 58 (7.14%) | 41 (4.24%) | 52 (5.78%) | 0.029 | 39 (3.17%) | 29 (2.56%) | 33 (2.57%) | 0.576 |
| AMI, stroke, or CV Death | 8.48% | 4.42% | 5.71% | 0.008 | 3.84% | 2.66% | 2.16% | 0.269 |
#Tables can be exported to CSV, HTML, LaTeX, PDF, Markdown, Word or Excel
#export2csv(restab, file='table1.csv'), exports to CSV format
#export2html(restab, file='table1.html'), exports to HTML format
#export2latex(restab, file='table1.tex'), exports to LaTeX format (to be included in Swaeave documents R chunks)
#export2pdf(restab, file='table1.pdf'), exports to PDF format
#export2md(restab, file='table1.md'), to be included inside Markdown documents R chunks
#export2word(restab, file='table1.docx'), exports to Word format
#export2xls(restab, file='table1.xlsx'), exports to Excel format
A new function has been implemented in the compareGroups package, which is called missingTable. This function applies to both compareGroups and createTable class objects. This last option is useful when the table is already created. To illustrate it, we will use the REGICOR data set, comparing missing rates of all variables by year:
# from a compareGroups object
data(regicor)
res <- compareGroups(year ~ . - id, regicor)
missingTable(res)
## Warning in cor(as.integer(x), as.integer(y)): the standard deviation is
## zero
## Warning in cor(as.integer(x), as.integer(y)): the standard deviation is
## zero
##
## --------Missingness table by 'year'---------
##
## ____________________________________________________
## 1995 2000 2005 p.overall
## N=431 N=786 N=1077
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## age 0 (0.00%) 0 (0.00%) 0 (0.00%) .
## sex 0 (0.00%) 0 (0.00%) 0 (0.00%) .
## smoker 16 (3.71%) 28 (3.56%) 17 (1.58%) 0.010
## sbp 3 (0.70%) 11 (1.40%) 0 (0.00%) <0.001
## dbp 3 (0.70%) 11 (1.40%) 0 (0.00%) <0.001
## histhtn 0 (0.00%) 0 (0.00%) 8 (0.74%) 0.015
## txhtn 0 (0.00%) 0 (0.00%) 43 (3.99%) <0.001
## chol 28 (6.50%) 71 (9.03%) 2 (0.19%) <0.001
## hdl 30 (6.96%) 38 (4.83%) 1 (0.09%) <0.001
## triglyc 28 (6.50%) 34 (4.33%) 1 (0.09%) <0.001
## ldl 43 (9.98%) 98 (12.5%) 27 (2.51%) <0.001
## histchol 0 (0.00%) 15 (1.91%) 6 (0.56%) 0.001
## txchol 0 (0.00%) 13 (1.65%) 42 (3.90%) <0.001
## height 8 (1.86%) 15 (1.91%) 12 (1.11%) 0.318
## weight 8 (1.86%) 15 (1.91%) 12 (1.11%) 0.318
## bmi 8 (1.86%) 15 (1.91%) 12 (1.11%) 0.318
## phyact 64 (14.8%) 22 (2.80%) 2 (0.19%) <0.001
## pcs 34 (7.89%) 123 (15.6%) 83 (7.71%) <0.001
## mcs 34 (7.89%) 123 (15.6%) 83 (7.71%) <0.001
## cv 33 (7.66%) 45 (5.73%) 53 (4.92%) 0.118
## tocv 33 (7.66%) 45 (5.73%) 53 (4.92%) 0.118
## death 44 (10.2%) 48 (6.11%) 54 (5.01%) 0.001
## todeath 44 (10.2%) 48 (6.11%) 54 (5.01%) 0.001
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
#summarize missingness
# the funcation counts the # of missing value for each variable
propmiss <- function(dataframe) {
m <- sapply(dataframe, function(x) {
data.frame(
nmiss=sum(is.na(x)),
n=length(x),
propmiss=sum(is.na(x))/length(x)
)
})
d <- data.frame(t(m))
d <- sapply(d, unlist)
d <- as.data.frame(d)
d$variable <- row.names(d)
row.names(d) <- NULL
d <- cbind(d[ncol(d)],d[-ncol(d)])
return(d[order(d$propmiss), ])
}
#replace NA to mean and mode
Mode<-function(x,na.rm){
xtab <- table(x)
xmode <- names(which(xtab == max(xtab)))
if (length(xmode) > 1) xmode <- ">1 mode"
return(xmode)
}
df_test
for (var in 1:ncol(df_test)) {
if (class(df_test[,var])=="numeric") {
df_test[is.na(df_test[,var]),var] <- mean(df_test[,var], na.rm = TRUE)
} else if (class(df_test[,var]) %in% c("character", "factor")) {
df_test[is.na(df_test[,var]),var] <- Mode(df_test[,var], na.rm = TRUE)
}
}
ac<-df_test
#load package with the data
#install.packages("devtools") # if not already installed
# library(devtools)
# install_git("https://github.com/ccolonescu/PoEdata")
#load data
data("cps4_small", package="PoEdata")
names(cps4_small)
## [1] "wage" "educ" "exper" "hrswk" "married" "female" "metro"
## [8] "midwest" "south" "west" "black" "asian"
#run the model
dnosouth <- cps4_small[which(cps4_small$south==0),]#no south
dsouth <- cps4_small[which(cps4_small$south==1),] #south
mod5ns <- lm(wage~educ+black*female, data=dnosouth)
mod5s <- lm(wage~educ+black*female, data=dsouth)
mod6 <- lm(wage~educ+black*female+south/(educ+black*female),
data=cps4_small)
#Stargazer cheat sheets
#https://www.jakeruss.com/cheatsheets/stargazer/
print(stargazer(mod6, mod5ns, mod5s, header=FALSE,
type="text",
title="Model comparison, 'wage' equation",
keep.stat="n",digits=2, single.row=TRUE,
intercept.bottom=FALSE))
##
## Model comparison, 'wage' equation
## ==================================================================
## Dependent variable:
## -----------------------------------------------
## wage
## (1) (2) (3)
## ------------------------------------------------------------------
## Constant -6.61*** (2.34) -6.61*** (2.30) -2.66 (3.42)
## educ 2.17*** (0.17) 2.17*** (0.16) 1.86*** (0.24)
## black -5.09* (2.64) -5.09* (2.60) -3.38 (2.58)
## female -5.01*** (0.90) -5.01*** (0.89) -4.10*** (1.58)
## south 3.94 (4.05)
## black:female 5.31 (3.50) 5.31 (3.45) 2.37 (3.38)
## educ:south -0.31 (0.29)
## black:south 1.70 (3.63)
## female:south 0.90 (1.77)
## black:female:south -2.94 (4.79)
## ------------------------------------------------------------------
## Observations 1,000 704 296
## ==================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
## [1] ""
## [2] "Model comparison, 'wage' equation"
## [3] "=================================================================="
## [4] " Dependent variable: "
## [5] " -----------------------------------------------"
## [6] " wage "
## [7] " (1) (2) (3) "
## [8] "------------------------------------------------------------------"
## [9] "Constant -6.61*** (2.34) -6.61*** (2.30) -2.66 (3.42) "
## [10] "educ 2.17*** (0.17) 2.17*** (0.16) 1.86*** (0.24) "
## [11] "black -5.09* (2.64) -5.09* (2.60) -3.38 (2.58) "
## [12] "female -5.01*** (0.90) -5.01*** (0.89) -4.10*** (1.58)"
## [13] "south 3.94 (4.05) "
## [14] "black:female 5.31 (3.50) 5.31 (3.45) 2.37 (3.38) "
## [15] "educ:south -0.31 (0.29) "
## [16] "black:south 1.70 (3.63) "
## [17] "female:south 0.90 (1.77) "
## [18] "black:female:south -2.94 (4.79) "
## [19] "------------------------------------------------------------------"
## [20] "Observations 1,000 704 296 "
## [21] "=================================================================="
## [22] "Note: *p<0.1; **p<0.05; ***p<0.01"
#output the results to html, we create a HTML file that can be easily copied/pasted into Word/Excel.
#HTML output may look like meaningless code in the console. use results = F to mask them.
# add "'" in front of the cell contents in excel to avoid function
xxx<-stargazer(mod6, mod5ns, mod5s, header=FALSE,
type="html",
title="Model comparison, 'wage' equation",
keep.stat="n",digits=2, single.row=TRUE,
intercept.bottom=FALSE, out= "xxx.html")
library(rio)
#export(as.data.frame(xxx), "xxx.csv")
#library(devtools)
#install_github("johnjosephhorton/JJHmisc")
library(JJHmisc)
library(help = JJHmisc)
lsf.str("package:JJHmisc")
## AddTableNote : function (stargazer.object, out.file, note, adjust = 3, table.type = "table")
## big.letter.theme : function (font.size)
## clusterSE : function (cluster.name, model, df)
## Codebook.Entry : function (var.name, label, note)
## createLaTeXparameters : function (key.list, parameter.file.location, append = FALSE)
## dzip : function (keys, values)
## fix.sample.size.line : function (table.line)
## flattenList : function (first.key, values)
## genParamAdder : function (parameter.file)
## get_field_names : function (con, tablename)
## get.line : function (linenumber, file)
## GetData : function (table.name, file.name, refresh)
## HasChange : function (x, k)
## insert.note : function (linenumber, note, file)
## kill.lines : function (start, end, file)
## latex.column.label : function (width, number, label, include.DV = TRUE)
## latex.row.label : function (width, label, rule.pts = 0, new.lines = 0, new.line.char = FALSE)
## NoteWrapper : function (x, text.width = 0.9)
## pp.big.number : function (e)
## regression.table : function (models, renames, out.file)
## render : function (x, lst)
## spanning.latex.column.label : function (num.columns, width, label)
## starLabel : function (p)
## tableContinuous2 : function (vars, weights = NA, subset = NA, group = NA, stats = c("n",
## "min", "q1", "median", "mean", "q3", "max", "s", "iqr", "na"),
## prec = 1, col.tit = NA, col.tit.font = c("bf", "", "sf", "it",
## "rm"), print.pval = c("none", "anova", "kruskal"), pval.bound = 10^-4,
## declare.zero = 10^-10, cap = "", lab = "", font.size = "footnotesize",
## longtable = TRUE, disp.cols = NA, nams = NA)
## tableNominal2 : function (vars, weights = NA, subset = NA, group = NA, miss.cat = NA,
## print.pval = c("none", "fisher", "chi2"), pval.bound = 10^-4,
## fisher.B = 2000, vertical = TRUE, cap = "", lab = "", col.tit.font = c("bf",
## "", "sf", "it", "rm"), font.size = "footnotesize", longtable = TRUE,
## nams = NA, cumsum = TRUE)
## wasItUsed : function (r.file, data.frame)
## write.image : function (filename, g, width = 5, height = 5, format = "png")
## writeImage : function (g, file.name.no.ext, width = 15, height = 8, path = "../../writeup/plots/",
## position = "h", line.width = "0.8", caption = "", notes = "",
## include.tex.wrapper = FALSE)
#read the source code of AddTableNote function
JJHmisc:::AddTableNote
## function (stargazer.object, out.file, note, adjust = 3, table.type = "table")
## {
## "Adds a table note to a stargazer table. It assumes that the\n last three lines of stargazer output are useless."
## n <- length(stargazer.object)
## write(stargazer.object[1:(n - adjust)], out.file)
## write(c("\\end{tabular}"), out.file, append = TRUE)
## write(note, out.file, append = TRUE)
## write(c(paste0("\\end{", table.type, "}")), out.file, append = TRUE)
## }
## <bytecode: 0x7fef18b20928>
## <environment: namespace:JJHmisc>
knitr::include_graphics("rFunctionsList.pdf")
Useful functions
library(rstan)
library(bayesplot)
stancode <- 'data {real y_mean;}
parameters {real y;}
model {y ~ normal(y_mean,1);}'
mod <- stan_model(model_code = stancode)
simple.fit <- sampling(mod, data = list(y_mean = 0))
##
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 7e-06 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.006417 seconds (Warm-up)
## Chain 1: 0.006959 seconds (Sampling)
## Chain 1: 0.013376 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 2e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.006231 seconds (Warm-up)
## Chain 2: 0.00617 seconds (Sampling)
## Chain 2: 0.012401 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 2e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.006364 seconds (Warm-up)
## Chain 3: 0.006119 seconds (Sampling)
## Chain 3: 0.012483 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 3e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.0061 seconds (Warm-up)
## Chain 4: 0.005758 seconds (Sampling)
## Chain 4: 0.011858 seconds (Total)
## Chain 4:
simple.fit2 <- sampling(mod, data = list(y_mean = 5))
##
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 3e-06 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.00611 seconds (Warm-up)
## Chain 1: 0.005841 seconds (Sampling)
## Chain 1: 0.011951 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 3e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.006205 seconds (Warm-up)
## Chain 2: 0.006822 seconds (Sampling)
## Chain 2: 0.013027 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 2e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.006055 seconds (Warm-up)
## Chain 3: 0.005907 seconds (Sampling)
## Chain 3: 0.011962 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 2e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.006037 seconds (Warm-up)
## Chain 4: 0.006337 seconds (Sampling)
## Chain 4: 0.012374 seconds (Total)
## Chain 4:
Comments are indicated by // in Stan. The write("model code", "file_name") bit allows us to write the Stan model in our R script and output the file to the working directory (or you can set a different file path).
write("// Stan model for simple linear regression
data {
int < lower = 1 > N; // Sample size
vector[N] x; // Predictor
vector[N] y; // Outcome
}
parameters {
real alpha; // Intercept
real beta; // Slope (regression coefficients)
real < lower = 0 > sigma; // Error SD
}
model {
y ~ normal(alpha + x * beta , sigma);
}
generated quantities {
} // The posterior predictive distribution",
"stan_model1.stan")
# stanc("stan_model1.stan")
stan_model1 <- "stan_model1.stan"
#load data
# Adding stringsAsFactors = F means that numeric variables won't be
# read in as factors/categorical variables
seaice <- read.csv("seaice.csv", stringsAsFactors = F)
head(seaice)
## year extent_north extent_south
## 1 1979 12.328 11.700
## 2 1980 12.337 11.230
## 3 1981 12.127 11.435
## 4 1982 12.447 11.640
## 5 1983 12.332 11.389
## 6 1984 11.910 11.454
#To explore the answer to that question, first we can make a figure.
plot(extent_north ~ year, pch = 20, data = seaice)
#preparing data for stan
x <- I(seaice$year - 1978) #rename the variables and index the years from 1 to 39.
y <- seaice$extent_north
N <- length(seaice$year)
stan_data <- list(N = N, x = x, y = y)
#fit the model with stan
stan_model1.fit <- stan(file = stan_model1,
data = stan_data,
warmup = 500, iter = 1000,
chains = 4, cores = 2, thin = 1)
stan_model1.fit
## Inference for Stan model: stan_model1.
## 4 chains, each with iter=1000; warmup=500; thin=1;
## post-warmup draws per chain=500, total post-warmup draws=2000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## alpha 12.55 0.00 0.08 12.40 12.50 12.55 12.61 12.69 849 1
## beta -0.05 0.00 0.00 -0.06 -0.06 -0.05 -0.05 -0.05 958 1
## sigma 0.23 0.00 0.03 0.18 0.21 0.22 0.24 0.28 904 1
## lp__ 37.48 0.04 1.19 34.32 36.91 37.80 38.36 38.87 797 1
##
## Samples were drawn using NUTS(diag_e) at Mon Jan 7 16:27:03 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
#look at the full posterior of our parameters by extracting them from the model object.
#extract() puts the posterior estimates for each parameter into a list.
posterior <- extract(stan_model1.fit)
str(posterior)
## List of 4
## $ alpha: num [1:2000(1d)] 12.5 12.5 12.5 12.4 12.6 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
## $ beta : num [1:2000(1d)] -0.0525 -0.0523 -0.0532 -0.05 -0.0563 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
## $ sigma: num [1:2000(1d)] 0.245 0.231 0.246 0.215 0.275 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
## $ lp__ : num [1:2000(1d)] 38.2 38.6 38.2 37.5 37 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
#We can also look at the posterior densities & histograms.
stan_dens(stan_model1.fit)
stan_hist(stan_model1.fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#we can generate plots which indicate the mean parameter estimates and
#any credible intervals we may be interested in.
#Note that the 95% credible intervals for the beta and
#sigma parameters are very small, thus you only see the dots.
plot(stan_model1.fit, show_density = FALSE, ci_level = 0.5,
outer_level = 0.95, fill_color = "salmon")
## ci_level: 0.5 (50% intervals)
## outer_level: 0.95 (95% intervals)
#comparing with lm function
lm1 <- lm(y ~ x)
summary(lm1)
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.49925 -0.17713 0.04898 0.16923 0.32829
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.555888 0.071985 174.4 <0.0000000000000002 ***
## x -0.054574 0.003137 -17.4 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2205 on 37 degrees of freedom
## Multiple R-squared: 0.8911, Adjusted R-squared: 0.8881
## F-statistic: 302.7 on 1 and 37 DF, p-value: < 0.00000000000000022
lm_alpha <- summary(lm1)$coeff[1] # the intercept
lm_beta <- summary(lm1)$coeff[2] # the slope
lm_sigma <- sigma(lm1) # the residual error
plot(y ~ x, pch = 20)
abline(lm1, col = 2, lty = 2, lw = 3)
abline( mean(posterior$alpha), mean(posterior$beta), col = 6, lw = 2)
The result is identical to the lm output. This is because we are using a simple model, and have put non-informative priors on our parameters.
One way to visualize the variability in our estimation of the regression line is to plot multiple estimates from the posterior.
library(gdata)
plot(y ~ x, pch = 20)
for (i in 1:500) {
abline(posterior$alpha[i], posterior$beta[i], col = "gray", lty = 1)
}
abline(mean(posterior$alpha), mean(posterior$beta), col = 6, lw = 2)
####################################################################
Let’s try again, but now with more informative priors for the relationship between sea ice and time. We’re going to use normal priors with small standard deviations. If we were to use normal priors with very large standard deviations (say 1000, or 10,000), they would act very similarly to uniform priors.
write("// Stan model for simple linear regression
data {
int < lower = 1 > N; // Sample size
vector[N] x; // Predictor
vector[N] y; // Outcome
}
parameters {
real alpha; // Intercept
real beta; // Slope (regression coefficients)
real < lower = 0 > sigma; // Error SD
}
model {
alpha ~ normal(10, 0.1);
beta ~ normal(1, 0.1);
y ~ normal(alpha + x * beta , sigma);
}
generated quantities {}",
"stan_model2.stan")
stan_model2 <- "stan_model2.stan"
fit2 <- stan(stan_model2, data = stan_data, warmup = 500, iter = 1000, chains = 4, cores = 2, thin = 1)
posterior2 <- extract(fit2)
plot(y ~ x, pch = 20)
#abline(alpha, beta, col = 4, lty = 2, lw = 2)
abline(mean(posterior2$alpha), mean(posterior2$beta), col = 3, lw = 2)
abline(mean(posterior$alpha), mean(posterior$beta), col = 36, lw = 3)
This is a common issue in Bayesian modelling, if your prior distributions are very narrow and yet don’t fit your understanding of the system or the distribution of your data, you could run models that do not meaningfully explain variation in your data. However, that isn’t to say that you shouldn’t choose somewhat informative priors, you do want to use previous analyses and understanding of your study system inform your model priors and design. You just need to think carefully about each modelling decision you make!
Convergence Diagnostics For traceplots, we can view them directly from the posterior:
plot(posterior$alpha, type = "l")
plot(posterior$beta, type = "l")
plot(posterior$sigma, type = "l")
traceplot(stan_model1.fit)
Example of poor convergence
Try running a model for only 50 iterations and check the traceplots.
fit_bad <- stan(stan_model1, data = stan_data, warmup = 25, iter = 50, chains = 4, cores = 2, thin = 1)
## Warning: There were 5 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low
## Warning: Examine the pairs() plot to diagnose sampling problems
posterior_bad <- extract(fit_bad)
plot(posterior_bad$alpha, type = "l")
plot(posterior_bad$beta, type = "l")
plot(posterior_bad$sigma, type = "l")
Parameter summaries
We can also get summaries of the parameters through the posterior directly. Let’s also plot the non-Bayesian linear model values to make sure our model is doing what we think it is…
par(mfrow = c(1,3))
plot(density(posterior$alpha), main = "Alpha")
abline(v = lm_alpha, col = 4, lty = 2)
plot(density(posterior$beta), main = "Beta")
abline(v = lm_beta, col = 4, lty = 2)
plot(density(posterior$sigma), main = "Sigma")
abline(v = lm_sigma, col = 4, lty = 2)
From the posterior we can directly calculate the probability of any parameter being over or under a certain value of interest.
#Probablility that beta is >0:
sum(posterior$beta>0)/length(posterior$beta)
## [1] 0
#Probablility that beta is >0.2:
sum(posterior$beta>0.2)/length(posterior$beta)
## [1] 0
Posterior Predictive Checks
For prediction and as another form of model diagnostic, Stan can use random number generators to generate predicted values for each data point, at each iteration. This way we can generate predictions that also represent the uncertainties in our model and our data generation process. We generate these using the Generated Quantities block. This block can be used to get any other information we want about the posterior, or make predictions for new data.
Note that vectorization is not supported in the GQ (generated quantities) block, so we have to put it in a loop. But since this is compiled to C++, loops are actually quite fast and Stan only evaluates the GQ block once per iteration, so it won’t add too much time to your sampling. Typically, the data generating functions will be the distributions you used in the model block but with an _rng suffix. (Double-check in the Stan manual to see which sampling statements have corresponding rng functions already coded up.)
write("// Stan model for simple linear regression
data {
int < lower = 1 > N; // Sample size
vector[N] x; // Predictor
vector[N] y; // Outcome
}
parameters {
real alpha; // Intercept
real beta; // Slope (regression coefficients)
real < lower = 0 > sigma; // Error SD
}
model {
y ~ normal(x * beta + alpha, sigma);
}
generated quantities {
real y_rep[N];
for (n in 1:N) {
y_rep[n] = normal_rng(x[n] * beta + alpha, sigma);
}
}",
"stan_model2_GQ.stan")
stan_model2_GQ <- "stan_model2_GQ.stan"
fit3 <- stan(stan_model2_GQ, data = stan_data, iter = 1000, chains = 4, cores = 2, thin = 1)
#Extracting the y_rep values from posterior.
#Each row is an iteration (single posterior estimate) from the model.
y_rep <- as.matrix(fit3, pars = "y_rep")
dim(y_rep)
## [1] 2000 39
#Comparing density of y with densities of y over 200 posterior draws.
bayesplot::ppc_dens_overlay(y, y_rep[1:200, ])
#We can also use this to compare estimates of summary statistics.
bayesplot::ppc_stat(y = y, yrep = y_rep, stat = "mean")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#We can investigate mean posterior prediction per datapoint vs
#the observed value for each datapoint (default line is 1:1)
bayesplot::ppc_scatter_avg(y = y, yrep = y_rep)
#Here is a list of currently available plots (bayesplot 1.2):
#available_ppc()
#can change the colour scheme in bayesplot too:
bayesplot::color_scheme_view(c("blue", "gray", "green", "pink", "purple",
"red","teal","yellow"))
#And you can even mix them:
color_scheme_view("mix-blue-red")
#You can set color schemes with:
color_scheme_set("blue")
data {
int<lower=0> N;
int<lower=0> T;
real x[T];
real y[N,T];
real xbar;
}
parameters {
real alpha[N];
real beta[N];
real mu_alpha;
real mu_beta; // beta.c in original bugs model
real<lower=0> sigmasq_y;
real<lower=0> sigmasq_alpha;
real<lower=0> sigmasq_beta;
}
transformed parameters {
real<lower=0> sigma_y; // sigma in original bugs model
real<lower=0> sigma_alpha;
real<lower=0> sigma_beta;
sigma_y = sqrt(sigmasq_y);
sigma_alpha = sqrt(sigmasq_alpha);
sigma_beta = sqrt(sigmasq_beta);
}
model {
mu_alpha ~ normal(0, 100);
mu_beta ~ normal(0, 100);
sigmasq_y ~ inv_gamma(0.001, 0.001);
sigmasq_alpha ~ inv_gamma(0.001, 0.001);
sigmasq_beta ~ inv_gamma(0.001, 0.001);
alpha ~ normal(mu_alpha, sigma_alpha); // vectorized
beta ~ normal(mu_beta, sigma_beta); // vectorized
for (n in 1:N)
for (t in 1:T)
y[n,t] ~ normal(alpha[n] + beta[n] * (x[t] - xbar), sigma_y);
}
generated quantities {
real alpha0;
alpha0 = mu_alpha - xbar * mu_beta;
}
library(tidyverse)
df <- read_delim("rats.txt", delim = " ")
## Parsed with column specification:
## cols(
## day8 = col_double(),
## day15 = col_double(),
## day22 = col_double(),
## day29 = col_double(),
## day36 = col_double()
## )
df
## # A tibble: 30 x 5
## day8 day15 day22 day29 day36
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 151 199 246 283 320
## 2 145 199 249 293 354
## 3 147 214 263 312 328
## 4 155 200 237 272 297
## 5 135 188 230 280 323
## 6 159 210 252 298 331
## 7 141 189 231 275 305
## 8 159 201 248 297 338
## 9 177 236 285 350 376
## 10 134 182 220 260 296
## # … with 20 more rows
y <- as.matrix(df)
y
## day8 day15 day22 day29 day36
## [1,] 151 199 246 283 320
## [2,] 145 199 249 293 354
## [3,] 147 214 263 312 328
## [4,] 155 200 237 272 297
## [5,] 135 188 230 280 323
## [6,] 159 210 252 298 331
## [7,] 141 189 231 275 305
## [8,] 159 201 248 297 338
## [9,] 177 236 285 350 376
## [10,] 134 182 220 260 296
## [11,] 160 208 261 313 352
## [12,] 143 188 220 273 314
## [13,] 154 200 244 289 325
## [14,] 171 221 270 326 358
## [15,] 163 216 242 281 312
## [16,] 160 207 248 288 324
## [17,] 142 187 234 280 316
## [18,] 156 203 243 283 317
## [19,] 157 212 259 307 336
## [20,] 152 203 246 286 321
## [21,] 154 205 253 298 334
## [22,] 139 190 225 267 302
## [23,] 146 191 229 272 302
## [24,] 157 211 250 285 323
## [25,] 132 185 237 286 331
## [26,] 160 207 257 303 345
## [27,] 169 216 261 295 333
## [28,] 157 205 248 289 316
## [29,] 137 180 219 258 291
## [30,] 153 200 244 286 324
x <- c(8,15,22,29,36)
xbar <- mean(x)
N <- nrow(y)
T <- ncol(y)
library(rstan)
options(mc.cores = parallel::detectCores())
#sample from that posterior distribution
rats_fit <- rstan::sampling(rats,
data = list(y,x,xbar,N,T))
rstan::summary(rats_fit)
## $summary
## mean se_mean sd 2.5%
## alpha[1] 239.9101176 0.032107296 2.6578237 234.5488751
## alpha[2] 247.7710919 0.029412775 2.6246964 242.6236101
## alpha[3] 252.4522863 0.031460374 2.6109330 247.2652729
## alpha[4] 232.5657126 0.030989675 2.7144837 227.1594045
## alpha[5] 231.6309075 0.034440381 2.7045632 226.3269739
## alpha[6] 249.7481589 0.030889327 2.6572009 244.4669857
## alpha[7] 228.7758177 0.032963901 2.7019954 223.4477749
## alpha[8] 248.3437472 0.033694844 2.7709166 243.0794437
## alpha[9] 283.3338365 0.034574598 2.7318716 277.8795001
## alpha[10] 219.2564254 0.036961968 2.6462226 214.0145444
## alpha[11] 258.2312268 0.036763933 2.7278808 252.8447814
## alpha[12] 228.1444671 0.037685651 2.7347577 222.7936056
## alpha[13] 242.3802764 0.034283086 2.7192027 237.1123343
## alpha[14] 268.2839619 0.035878255 2.7584002 262.7762843
## alpha[15] 242.8368719 0.033333623 2.6940760 237.7157202
## alpha[16] 245.3301513 0.031967736 2.7005467 239.9721467
## alpha[17] 232.1730819 0.036344901 2.7550520 226.7701743
## alpha[18] 240.5326904 0.035837571 2.7692615 235.0787279
## alpha[19] 253.7550682 0.034482225 2.6786253 248.4006413
## alpha[20] 241.6498238 0.034534730 2.7206198 236.2303313
## alpha[21] 248.6097520 0.033301156 2.7067765 243.4666249
## alpha[22] 225.2224625 0.036639060 2.6407778 220.1035482
## alpha[23] 228.5369983 0.037710436 2.6966543 223.1865699
## alpha[24] 245.1108964 0.032851918 2.6535119 239.9802538
## alpha[25] 234.5014763 0.030387887 2.7142259 229.0450572
## alpha[26] 253.9817892 0.034379196 2.7662874 248.4171417
## alpha[27] 254.3184138 0.033148950 2.6496400 249.2526166
## alpha[28] 242.9907623 0.034766967 2.6758857 237.5824858
## alpha[29] 217.9332278 0.035070957 2.7041637 212.7331971
## alpha[30] 241.4534414 0.034657623 2.7229125 236.2199885
## beta[1] 6.0644501 0.003317410 0.2450428 5.5844573
## beta[2] 7.0528794 0.003605541 0.2557732 6.5299716
## beta[3] 6.4840340 0.003160061 0.2529363 5.9827008
## beta[4] 5.3421506 0.003280252 0.2528233 4.8292343
## beta[5] 6.5693428 0.003064760 0.2402844 6.0875833
## beta[6] 6.1741715 0.003276169 0.2432543 5.6955275
## beta[7] 5.9754878 0.003019117 0.2421691 5.4946117
## beta[8] 6.4204992 0.003119430 0.2468635 5.9338214
## beta[9] 7.0527735 0.003306984 0.2601278 6.5253891
## beta[10] 5.8549368 0.002974570 0.2434404 5.3741981
## beta[11] 6.7986044 0.003260923 0.2556692 6.2833233
## beta[12] 6.1200303 0.002905360 0.2418928 5.6448467
## beta[13] 6.1684915 0.003273787 0.2495856 5.6505664
## beta[14] 6.6907145 0.002901221 0.2580016 6.1858630
## beta[15] 5.4163624 0.003153302 0.2442228 4.9372132
## beta[16] 5.9217063 0.003038135 0.2404501 5.4428470
## beta[17] 6.2763483 0.003103078 0.2465650 5.7935908
## beta[18] 5.8473081 0.002921502 0.2426703 5.3645122
## beta[19] 6.4037304 0.003107552 0.2449235 5.9266222
## beta[20] 6.0547395 0.003046638 0.2448963 5.5626759
## beta[21] 6.4024657 0.003161328 0.2449800 5.9147997
## beta[22] 5.8627274 0.003171593 0.2435659 5.3924864
## beta[23] 5.7456716 0.003123035 0.2516111 5.2683213
## beta[24] 5.8917624 0.002771738 0.2472508 5.3973264
## beta[25] 6.9061627 0.003119167 0.2457827 6.4230444
## beta[26] 6.5481125 0.003068555 0.2427410 6.0677120
## beta[27] 5.9031704 0.002825252 0.2389190 5.4280174
## beta[28] 5.8487791 0.002891054 0.2440778 5.3605305
## beta[29] 5.6742810 0.003277911 0.2464774 5.1871659
## beta[30] 6.1342841 0.002726781 0.2374156 5.6732897
## mu_alpha 242.4769641 0.036647467 2.7496143 236.9904708
## mu_beta 6.1873138 0.001650494 0.1083687 5.9689339
## sigmasq_y 37.4509648 0.111887922 5.8444217 27.3319211
## sigmasq_alpha 217.8323501 1.135202289 64.6337673 124.8698817
## sigmasq_beta 0.2735714 0.001689988 0.1000803 0.1315568
## sigma_y 6.1015115 0.009067297 0.4717821 5.2279940
## sigma_alpha 14.6086201 0.034923241 2.1027775 11.1745193
## sigma_beta 0.5149685 0.001535668 0.0915475 0.3627077
## alpha0 106.3560610 0.052391111 3.6219112 99.0589096
## lp__ -438.6669450 0.216886901 7.0499385 -453.4130373
## 25% 50% 75% 97.5% n_eff
## alpha[1] 238.1675377 239.9333958 241.6984956 245.0355233 6852.434
## alpha[2] 246.0733276 247.7323203 249.5515871 253.0381621 7963.173
## alpha[3] 250.7041369 252.4513488 254.2281914 257.5878145 6887.533
## alpha[4] 230.7353875 232.5541032 234.3574918 237.8531382 7672.562
## alpha[5] 229.7891779 231.5972158 233.4234115 236.8976075 6166.779
## alpha[6] 247.9537566 249.7786087 251.5766478 254.8018159 7400.003
## alpha[7] 226.9337347 228.7873045 230.6294379 234.0227400 6718.804
## alpha[8] 246.4098075 248.2996921 250.2499589 253.9513876 6762.698
## alpha[9] 281.5543810 283.3615450 285.2205777 288.5110862 6243.186
## alpha[10] 217.4248965 219.2374032 221.0676174 224.4912249 5125.575
## alpha[11] 256.4146007 258.2244878 260.0467161 263.6088276 5505.628
## alpha[12] 226.3276403 228.2052378 229.9760624 233.5149403 5266.058
## alpha[13] 240.5610480 242.3535576 244.2290557 247.6511055 6291.053
## alpha[14] 266.4408861 268.3543095 270.1169785 273.6909386 5910.877
## alpha[15] 241.0688565 242.8195377 244.6080161 248.1623050 6532.127
## alpha[16] 243.5343836 245.3560795 247.1501472 250.4967635 7136.407
## alpha[17] 230.2873536 232.2097496 234.0473887 237.4544068 5746.093
## alpha[18] 238.7007733 240.5132073 242.3822496 246.0491638 5971.051
## alpha[19] 251.9807564 253.7744143 255.5407144 259.0320223 6034.390
## alpha[20] 239.8380285 241.6809851 243.4791093 247.0318862 6206.169
## alpha[21] 246.7655637 248.5731640 250.4145555 254.1298452 6606.724
## alpha[22] 223.4544631 225.2162203 226.9496564 230.3653646 5194.875
## alpha[23] 226.7496189 228.5775870 230.3223556 233.7823100 5113.608
## alpha[24] 243.3497722 245.1177653 246.8428643 250.3437867 6524.100
## alpha[25] 232.7671886 234.5388923 236.2405042 239.7228800 7977.944
## alpha[26] 252.1372000 253.9935830 255.8300567 259.4362636 6474.454
## alpha[27] 252.4940461 254.3519698 256.1306437 259.5492879 6389.019
## alpha[28] 241.1898075 242.9504615 244.7909839 248.3932803 5923.815
## alpha[29] 216.1015136 217.8699391 219.7309097 223.3787407 5945.259
## alpha[30] 239.6092495 241.4320058 243.2925863 246.7302783 6172.624
## beta[1] 5.9016868 6.0602342 6.2295013 6.5420962 5456.140
## beta[2] 6.8826410 7.0528893 7.2319203 7.5503089 5032.332
## beta[3] 6.3180118 6.4857643 6.6544347 6.9711355 6406.654
## beta[4] 5.1757242 5.3423518 5.5111674 5.8414451 5940.455
## beta[5] 6.4076781 6.5706441 6.7288095 7.0411539 6146.929
## beta[6] 6.0091246 6.1776375 6.3427587 6.6495631 5513.005
## beta[7] 5.8158684 5.9753670 6.1360390 6.4496104 6433.946
## beta[8] 6.2573140 6.4210072 6.5842162 6.9064715 6262.721
## beta[9] 6.8829068 7.0492251 7.2239928 7.5697179 6187.417
## beta[10] 5.6963592 5.8553535 6.0152041 6.3303438 6697.871
## beta[11] 6.6320326 6.8038259 6.9667485 7.3088905 6147.175
## beta[12] 5.9624477 6.1203819 6.2827176 6.5943635 6931.800
## beta[13] 6.0050765 6.1685567 6.3315431 6.6565623 5812.169
## beta[14] 6.5156479 6.6855874 6.8674784 7.1907583 7908.303
## beta[15] 5.2558817 5.4176733 5.5789920 5.9111793 5998.478
## beta[16] 5.7599369 5.9215235 6.0876286 6.3764558 6263.768
## beta[17] 6.1162209 6.2763330 6.4426962 6.7593673 6313.605
## beta[18] 5.6871604 5.8473301 6.0057197 6.3329377 6899.549
## beta[19] 6.2328684 6.4047654 6.5708075 6.8761905 6211.893
## beta[20] 5.8913325 6.0545078 6.2197273 6.5378323 6461.342
## beta[21] 6.2385311 6.4065739 6.5686046 6.8716331 6005.127
## beta[22] 5.6945705 5.8644918 6.0289380 6.3420863 5897.635
## beta[23] 5.5691256 5.7417880 5.9176685 6.2343153 6490.916
## beta[24] 5.7303897 5.8956260 6.0487064 6.3805677 7957.397
## beta[25] 6.7425226 6.9051803 7.0704856 7.3998105 6209.051
## beta[26] 6.3902110 6.5469365 6.7130224 7.0243083 6257.755
## beta[27] 5.7390866 5.9034882 6.0647555 6.3668808 7151.333
## beta[28] 5.6875684 5.8456505 6.0137796 6.3366710 7127.613
## beta[29] 5.5104588 5.6764913 5.8348260 6.1467949 5654.051
## beta[30] 5.9759472 6.1376583 6.2936226 6.6079329 7580.853
## mu_alpha 240.6768502 242.5065350 244.3402045 247.6967971 5629.316
## mu_beta 6.1156840 6.1904349 6.2572214 6.4017220 4311.020
## sigmasq_y 33.2734113 37.0175766 41.0758986 50.2616686 2728.454
## sigmasq_alpha 171.5328956 206.8671150 251.6332633 376.1496098 3241.697
## sigmasq_beta 0.2021241 0.2569905 0.3234950 0.5170319 3506.953
## sigma_y 5.7683110 6.0842071 6.4090482 7.0895462 2707.243
## sigma_alpha 13.0970568 14.3828757 15.8629525 19.3945767 3625.414
## sigma_beta 0.4495821 0.5069423 0.5687662 0.7190493 3553.844
## alpha0 103.9324694 106.3789046 108.7889424 113.3425357 4779.257
## lp__ -443.1543728 -438.4595874 -433.7743000 -425.7954482 1056.584
## Rhat
## alpha[1] 1.0000846
## alpha[2] 0.9992497
## alpha[3] 0.9996211
## alpha[4] 0.9996334
## alpha[5] 0.9992586
## alpha[6] 0.9996418
## alpha[7] 0.9993552
## alpha[8] 1.0001294
## alpha[9] 0.9994213
## alpha[10] 1.0003155
## alpha[11] 0.9996596
## alpha[12] 0.9995462
## alpha[13] 0.9992462
## alpha[14] 0.9995434
## alpha[15] 0.9996318
## alpha[16] 0.9991470
## alpha[17] 0.9993781
## alpha[18] 0.9995428
## alpha[19] 0.9991064
## alpha[20] 0.9998488
## alpha[21] 1.0001483
## alpha[22] 1.0000499
## alpha[23] 0.9999544
## alpha[24] 0.9993894
## alpha[25] 1.0000089
## alpha[26] 0.9994473
## alpha[27] 0.9999175
## alpha[28] 0.9991887
## alpha[29] 0.9996949
## alpha[30] 0.9993756
## beta[1] 0.9997811
## beta[2] 0.9993809
## beta[3] 0.9992219
## beta[4] 0.9995006
## beta[5] 0.9995295
## beta[6] 1.0000908
## beta[7] 0.9995685
## beta[8] 0.9994455
## beta[9] 0.9992980
## beta[10] 0.9995895
## beta[11] 0.9994011
## beta[12] 0.9994513
## beta[13] 1.0001090
## beta[14] 0.9997021
## beta[15] 0.9995223
## beta[16] 0.9993061
## beta[17] 0.9994898
## beta[18] 0.9993602
## beta[19] 1.0001532
## beta[20] 0.9997914
## beta[21] 0.9995996
## beta[22] 0.9999703
## beta[23] 0.9992240
## beta[24] 0.9991987
## beta[25] 0.9996382
## beta[26] 0.9995108
## beta[27] 0.9998728
## beta[28] 0.9998624
## beta[29] 0.9993781
## beta[30] 0.9992577
## mu_alpha 0.9996557
## mu_beta 0.9994931
## sigmasq_y 1.0005936
## sigmasq_alpha 1.0015092
## sigmasq_beta 1.0001771
## sigma_y 1.0005634
## sigma_alpha 1.0012839
## sigma_beta 1.0001317
## alpha0 0.9996848
## lp__ 1.0021535
##
## $c_summary
## , , chains = chain:1
##
## stats
## parameter mean sd 2.5% 25%
## alpha[1] 239.7972483 2.60639912 234.7125902 238.0280946
## alpha[2] 247.8251074 2.64178150 242.6955335 246.0865869
## alpha[3] 252.4075042 2.42149759 247.8130702 250.6569773
## alpha[4] 232.5820993 2.78701190 227.1576092 230.7455688
## alpha[5] 231.6376926 2.84620844 226.1038192 229.6695463
## alpha[6] 249.8192304 2.53488114 245.0315400 248.0620879
## alpha[7] 228.7994238 2.63838802 223.4539937 227.0364576
## alpha[8] 248.5014452 2.74863684 243.1739216 246.6413974
## alpha[9] 283.2586783 2.81259040 277.8128621 281.4368390
## alpha[10] 219.1674360 2.67444720 214.0436032 217.2310065
## alpha[11] 258.1701723 2.64754811 253.0265912 256.3253384
## alpha[12] 228.1554644 2.79487764 222.7650150 226.3199139
## alpha[13] 242.3428627 2.65700499 237.0236869 240.5946878
## alpha[14] 268.2743853 2.79573460 262.7344352 266.4646927
## alpha[15] 242.9084497 2.59779595 237.8332553 241.2410825
## alpha[16] 245.3464185 2.77616064 239.8875378 243.4896350
## alpha[17] 232.1369169 2.67040252 226.7696592 230.2940502
## alpha[18] 240.4450581 2.83306438 234.8916814 238.4722564
## alpha[19] 253.7630544 2.67810152 248.4073626 251.9712172
## alpha[20] 241.5329743 2.63624174 236.5329043 239.6132601
## alpha[21] 248.7009847 2.79779759 243.2702018 246.8070504
## alpha[22] 225.1551859 2.55138412 220.0214701 223.5254486
## alpha[23] 228.4582150 2.66217730 223.5799146 226.6417695
## alpha[24] 245.0845830 2.64331907 240.1111242 243.3219981
## alpha[25] 234.4808369 2.65146943 228.9687418 232.8232873
## alpha[26] 253.9711000 2.71312087 248.4631016 252.2982053
## alpha[27] 254.4157853 2.51900937 249.3283539 252.6990020
## alpha[28] 243.0326897 2.50076114 238.0271618 241.3584012
## alpha[29] 217.9333756 2.63575985 212.8333935 216.1945747
## alpha[30] 241.4935969 2.67247409 236.3254160 239.6412567
## beta[1] 6.0725365 0.25375758 5.5751724 5.9005951
## beta[2] 7.0572707 0.26576382 6.5272554 6.8768346
## beta[3] 6.4785965 0.25848103 5.9545301 6.3103362
## beta[4] 5.3378774 0.25591881 4.8546500 5.1632625
## beta[5] 6.5660464 0.23790468 6.1005313 6.4013382
## beta[6] 6.1681612 0.25476636 5.6588987 5.9948806
## beta[7] 5.9794942 0.23094359 5.5484815 5.8200887
## beta[8] 6.4250619 0.23607641 5.9714962 6.2638567
## beta[9] 7.0479970 0.26126136 6.5210483 6.8780313
## beta[10] 5.8554059 0.23998411 5.3694576 5.6999254
## beta[11] 6.8014372 0.25903763 6.2636761 6.6329094
## beta[12] 6.1278073 0.24319484 5.6680165 5.9540386
## beta[13] 6.1621114 0.23798081 5.6718732 6.0039676
## beta[14] 6.6924534 0.25517657 6.1751998 6.5243266
## beta[15] 5.4117894 0.25487922 4.9118082 5.2518834
## beta[16] 5.9153507 0.24062944 5.4173267 5.7554914
## beta[17] 6.2856552 0.23759743 5.8219599 6.1238942
## beta[18] 5.8436079 0.23358626 5.3756507 5.6839407
## beta[19] 6.3917618 0.24054249 5.9459660 6.2242411
## beta[20] 6.0612853 0.24378011 5.5561171 5.8987808
## beta[21] 6.4050348 0.23185595 5.9690453 6.2434648
## beta[22] 5.8692426 0.24153703 5.3892752 5.7122298
## beta[23] 5.7450932 0.25431074 5.2639867 5.5627092
## beta[24] 5.8953821 0.24383241 5.4019218 5.7432739
## beta[25] 6.9151386 0.24303343 6.4107486 6.7615367
## beta[26] 6.5509966 0.24114268 6.0755573 6.3911304
## beta[27] 5.9091451 0.23952422 5.4249514 5.7354093
## beta[28] 5.8409380 0.24879376 5.3467413 5.6777638
## beta[29] 5.6679922 0.25588099 5.1568610 5.4990600
## beta[30] 6.1339759 0.22812373 5.6922933 5.9931439
## mu_alpha 242.4409566 2.72575406 237.0569237 240.5731728
## mu_beta 6.1888002 0.10869891 5.9644545 6.1195778
## sigmasq_y 37.2086601 5.75626152 27.9376150 33.0239414
## sigmasq_alpha 220.0866857 65.32654435 126.9929368 172.8760351
## sigmasq_beta 0.2779605 0.10214142 0.1356473 0.2026671
## sigma_y 6.0820942 0.46584020 5.2856021 5.7466461
## sigma_alpha 14.6857275 2.10250193 11.2691142 13.1482331
## sigma_beta 0.5190288 0.09261822 0.3683032 0.4501857
## alpha0 106.2873529 3.59541168 99.1870589 103.8453352
## lp__ -438.2808538 6.93621723 -452.9071952 -442.7341882
## stats
## parameter 50% 75% 97.5%
## alpha[1] 239.8667947 241.5703332 244.7670055
## alpha[2] 247.7289901 249.5734675 253.1452331
## alpha[3] 252.3328409 254.0046945 257.0674777
## alpha[4] 232.5285156 234.3770633 238.0928150
## alpha[5] 231.7535331 233.5433699 237.1264280
## alpha[6] 249.8112129 251.5716195 254.5118904
## alpha[7] 228.8476008 230.6032508 233.8660793
## alpha[8] 248.3551000 250.3753860 254.0434512
## alpha[9] 283.2527414 285.1978995 288.6836865
## alpha[10] 219.1490321 221.1502726 224.2661844
## alpha[11] 258.1900953 259.9643330 263.4089542
## alpha[12] 228.1879610 229.9203585 233.6582897
## alpha[13] 242.2845373 244.1705987 247.5685440
## alpha[14] 268.3511802 270.0927333 273.8345933
## alpha[15] 242.9757276 244.4742942 248.2689828
## alpha[16] 245.4705045 247.2927467 250.5375328
## alpha[17] 232.1878184 234.0007558 237.1797169
## alpha[18] 240.4619626 242.4841965 245.6027199
## alpha[19] 253.7804844 255.6277787 259.0771619
## alpha[20] 241.5009428 243.3476970 246.8651783
## alpha[21] 248.6055741 250.5588138 254.3670005
## alpha[22] 225.1700510 226.8390154 229.8768042
## alpha[23] 228.3636389 230.2150977 233.7860752
## alpha[24] 245.0909185 246.8447215 250.0440502
## alpha[25] 234.5708732 236.1869283 239.5833601
## alpha[26] 253.9671944 255.7759811 259.4650977
## alpha[27] 254.4697205 256.1536545 259.1550804
## alpha[28] 243.0038281 244.7258425 248.0647944
## alpha[29] 217.8881224 219.5746914 223.2441091
## alpha[30] 241.4579842 243.2448924 246.5852838
## beta[1] 6.0752151 6.2493067 6.5823282
## beta[2] 7.0645709 7.2425927 7.5683913
## beta[3] 6.4803484 6.6487742 6.9810153
## beta[4] 5.3360573 5.5092392 5.8555246
## beta[5] 6.5677169 6.7122138 7.0615434
## beta[6] 6.1825373 6.3485085 6.6462179
## beta[7] 5.9774969 6.1335177 6.4464679
## beta[8] 6.4252280 6.5826515 6.9044311
## beta[9] 7.0467226 7.2209842 7.5427927
## beta[10] 5.8590136 6.0134038 6.3056013
## beta[11] 6.8035609 6.9767230 7.3106500
## beta[12] 6.1312030 6.2872747 6.5947122
## beta[13] 6.1648139 6.3202458 6.6258229
## beta[14] 6.6889807 6.8719880 7.1793983
## beta[15] 5.4098748 5.5856729 5.9164090
## beta[16] 5.9149716 6.0831396 6.3817143
## beta[17] 6.2803634 6.4473034 6.7484575
## beta[18] 5.8530618 5.9988373 6.2735124
## beta[19] 6.3896452 6.5594383 6.8673249
## beta[20] 6.0612798 6.2280326 6.5398203
## beta[21] 6.4044263 6.5531506 6.8705243
## beta[22] 5.8692256 6.0261199 6.3153091
## beta[23] 5.7492518 5.9173904 6.2195379
## beta[24] 5.8967558 6.0503862 6.3795245
## beta[25] 6.9226011 7.0701602 7.3823138
## beta[26] 6.5478469 6.7157076 7.0073336
## beta[27] 5.9081767 6.0759034 6.3783507
## beta[28] 5.8415502 6.0020651 6.3185916
## beta[29] 5.6672162 5.8308601 6.1617017
## beta[30] 6.1387495 6.2678595 6.5970509
## mu_alpha 242.4736702 244.3411848 247.6314297
## mu_beta 6.1911350 6.2561357 6.4063237
## sigmasq_y 36.7446865 40.7849134 49.9108488
## sigmasq_alpha 210.0478600 253.4239247 385.0172292
## sigmasq_beta 0.2639536 0.3309054 0.5193950
## sigma_y 6.0617396 6.3863067 7.0647606
## sigma_alpha 14.4930278 15.9192941 19.6218359
## sigma_beta 0.5137641 0.5752438 0.7206901
## alpha0 106.2827742 108.7140118 112.9597328
## lp__ -438.1609888 -433.2515897 -426.1289985
##
## , , chains = chain:2
##
## stats
## parameter mean sd 2.5% 25%
## alpha[1] 240.0747142 2.67480931 234.6014977 238.3510757
## alpha[2] 247.7791483 2.48076409 242.8598396 246.1585024
## alpha[3] 252.5400268 2.73344171 246.9316964 250.8053271
## alpha[4] 232.5043000 2.73911279 227.0104892 230.6651318
## alpha[5] 231.5724667 2.79956394 226.1547100 229.7264167
## alpha[6] 249.6979044 2.66348224 244.2892738 247.9633249
## alpha[7] 228.7417049 2.78991710 223.4282040 226.8823485
## alpha[8] 248.3748856 2.82742957 242.7569469 246.2869964
## alpha[9] 283.4186550 2.67413497 278.2628846 281.6506620
## alpha[10] 219.4607838 2.61574465 214.5820323 217.7577664
## alpha[11] 258.2763834 2.71785181 253.0133848 256.4410328
## alpha[12] 228.0904648 2.70402805 222.7266779 226.3245153
## alpha[13] 242.4535014 2.70159458 237.2739232 240.6666538
## alpha[14] 268.3856271 2.84523616 262.8407007 266.4889397
## alpha[15] 242.7291418 2.66002406 237.7145691 241.0109200
## alpha[16] 245.3076640 2.70145230 239.9451239 243.5293278
## alpha[17] 232.1204437 2.65399778 227.1262676 230.2473892
## alpha[18] 240.5630063 2.71429868 235.2813937 238.6986005
## alpha[19] 253.7622650 2.80493841 248.1296828 251.7788658
## alpha[20] 241.6292937 2.74220110 236.1444055 239.8148044
## alpha[21] 248.5653116 2.65089624 243.5092800 246.8583196
## alpha[22] 225.2495660 2.62349610 220.4590882 223.4004887
## alpha[23] 228.5922810 2.83136168 223.1915052 226.6264426
## alpha[24] 245.1382897 2.59282809 240.1975897 243.4616159
## alpha[25] 234.6212015 2.72561285 229.3240090 232.8437234
## alpha[26] 254.0276301 2.85061782 248.7552747 251.9345470
## alpha[27] 254.2933415 2.73820219 249.2976326 252.1929444
## alpha[28] 242.9818848 2.72103431 237.7555241 241.0687728
## alpha[29] 217.9564911 2.77248359 212.8217920 215.9531638
## alpha[30] 241.4062730 2.78460016 235.9520482 239.6014886
## beta[1] 6.0570757 0.23995987 5.5592126 5.9031738
## beta[2] 7.0455736 0.24823125 6.5486288 6.8808357
## beta[3] 6.4872947 0.24987883 5.9958932 6.3317133
## beta[4] 5.3374653 0.25112179 4.8184346 5.1709624
## beta[5] 6.5673969 0.23806505 6.0984647 6.4075377
## beta[6] 6.1782835 0.24368123 5.6953206 6.0084617
## beta[7] 5.9818976 0.24710363 5.4850908 5.8212416
## beta[8] 6.4203687 0.25324816 5.9238352 6.2552395
## beta[9] 7.0604632 0.25859440 6.5432910 6.8856921
## beta[10] 5.8493531 0.25098768 5.3367616 5.6847870
## beta[11] 6.7959306 0.24151337 6.3082028 6.6425787
## beta[12] 6.1163597 0.23988947 5.6477409 5.9607595
## beta[13] 6.1630874 0.24623088 5.6526444 6.0065985
## beta[14] 6.6914393 0.26509687 6.1824109 6.4991512
## beta[15] 5.4167923 0.23699726 4.9668301 5.2661803
## beta[16] 5.9246342 0.23004250 5.4726913 5.7701535
## beta[17] 6.2739346 0.24434763 5.7928751 6.1165333
## beta[18] 5.8491454 0.25495552 5.3555984 5.6817647
## beta[19] 6.4124240 0.24349639 5.9645970 6.2416357
## beta[20] 6.0550317 0.24775494 5.5418660 5.8954408
## beta[21] 6.3931212 0.24361230 5.9003630 6.2406859
## beta[22] 5.8615080 0.25233071 5.3859543 5.6697783
## beta[23] 5.7491811 0.24152150 5.2846098 5.5801701
## beta[24] 5.8915084 0.23950650 5.4132766 5.7417108
## beta[25] 6.9023390 0.24699475 6.4208513 6.7349543
## beta[26] 6.5462247 0.23979237 6.0720912 6.3912508
## beta[27] 5.8975327 0.24375454 5.4440496 5.7367358
## beta[28] 5.8544269 0.24488914 5.3743382 5.6923970
## beta[29] 5.6730817 0.23902089 5.2002456 5.5119920
## beta[30] 6.1310352 0.24644089 5.6468563 5.9620273
## mu_alpha 242.4516883 2.91564888 236.5638855 240.5695985
## mu_beta 6.1864800 0.10850239 5.9719642 6.1142867
## sigmasq_y 37.5397658 5.77165671 26.7517692 33.5337760
## sigmasq_alpha 221.7872129 70.26058358 123.6443747 171.6186410
## sigmasq_beta 0.2754960 0.09896423 0.1314987 0.2026667
## sigma_y 6.1090332 0.46872082 5.1722112 5.7908354
## sigma_alpha 14.7202996 2.25944628 11.1195484 13.1003297
## sigma_beta 0.5169419 0.09096892 0.3626275 0.4501852
## alpha0 106.3491284 3.69921527 98.8126777 103.9975807
## lp__ -438.8659929 7.16304424 -455.2337088 -443.4893529
## stats
## parameter 50% 75% 97.5%
## alpha[1] 240.1610520 241.7474582 245.3130745
## alpha[2] 247.7859673 249.4398650 252.9181520
## alpha[3] 252.6512448 254.4022858 257.8744096
## alpha[4] 232.4895592 234.2154667 238.0307830
## alpha[5] 231.5143614 233.4810615 237.3391664
## alpha[6] 249.7116914 251.4782609 254.6445344
## alpha[7] 228.7621856 230.6846605 234.2117680
## alpha[8] 248.4956397 250.4086693 253.9158397
## alpha[9] 283.3664683 285.3723334 288.5136279
## alpha[10] 219.3347147 221.2365367 224.7189544
## alpha[11] 258.2335374 260.1207969 263.4907740
## alpha[12] 228.2126934 229.9389770 233.2161719
## alpha[13] 242.4077600 244.2092457 247.7374655
## alpha[14] 268.3807680 270.2517317 273.8917225
## alpha[15] 242.6550634 244.5996129 247.8651141
## alpha[16] 245.3410139 247.0739714 250.3242025
## alpha[17] 232.0843693 233.9659784 237.2120838
## alpha[18] 240.5114448 242.4961966 246.0306898
## alpha[19] 253.8131664 255.5839424 259.2131635
## alpha[20] 241.6368721 243.5353991 247.2077761
## alpha[21] 248.5596940 250.3542836 253.4632007
## alpha[22] 225.1339228 227.0613126 230.3345192
## alpha[23] 228.7652007 230.5233440 233.7702358
## alpha[24] 245.1432921 246.7210075 250.3509015
## alpha[25] 234.6890684 236.4008436 239.7644741
## alpha[26] 253.9966584 256.0318602 259.1694514
## alpha[27] 254.2857658 256.2760987 259.6261167
## alpha[28] 242.9667833 244.8326251 248.6275956
## alpha[29] 217.8886642 219.9075268 223.7470454
## alpha[30] 241.3483505 243.3118613 246.5809743
## beta[1] 6.0574730 6.2188376 6.5209861
## beta[2] 7.0366741 7.2100928 7.5447620
## beta[3] 6.4893272 6.6537480 6.9565601
## beta[4] 5.3413451 5.5115431 5.8213470
## beta[5] 6.5620287 6.7306732 7.0246114
## beta[6] 6.1846659 6.3488964 6.6354444
## beta[7] 5.9792160 6.1429901 6.4750550
## beta[8] 6.4290286 6.5825511 6.8986123
## beta[9] 7.0561164 7.2312302 7.5773047
## beta[10] 5.8499079 6.0060804 6.3662178
## beta[11] 6.8022274 6.9474881 7.2997560
## beta[12] 6.1164113 6.2832467 6.5691386
## beta[13] 6.1658161 6.3144293 6.6566333
## beta[14] 6.6930047 6.8726290 7.1978639
## beta[15] 5.4104574 5.5675423 5.9054396
## beta[16] 5.9223063 6.0792712 6.3737321
## beta[17] 6.2783930 6.4368432 6.7621670
## beta[18] 5.8465627 6.0183761 6.3533501
## beta[19] 6.4079061 6.5856453 6.8731707
## beta[20] 6.0619902 6.2106457 6.5224253
## beta[21] 6.3977751 6.5564404 6.8495474
## beta[22] 5.8563616 6.0432677 6.3627721
## beta[23] 5.7498557 5.9144709 6.2191615
## beta[24] 5.8958496 6.0408314 6.3692476
## beta[25] 6.8971236 7.0724173 7.3810687
## beta[26] 6.5425162 6.7025793 7.0243543
## beta[27] 5.8968628 6.0601252 6.3730498
## beta[28] 5.8512035 6.0193571 6.3423237
## beta[29] 5.6722369 5.8355217 6.1295456
## beta[30] 6.1397007 6.2937324 6.6235798
## mu_alpha 242.4641938 244.4273035 247.8313812
## mu_beta 6.1887881 6.2585127 6.3987412
## sigmasq_y 37.1891048 41.1386489 50.0381626
## sigmasq_alpha 207.7786588 256.6538142 398.3097495
## sigmasq_beta 0.2586505 0.3282218 0.4925724
## sigma_y 6.0982870 6.4139417 7.0737654
## sigma_alpha 14.4145293 16.0204179 19.9576989
## sigma_beta 0.5085770 0.5729063 0.7018350
## alpha0 106.3449729 108.8658653 113.2457998
## lp__ -438.4862993 -433.9871375 -425.5198251
##
## , , chains = chain:3
##
## stats
## parameter mean sd 2.5% 25%
## alpha[1] 239.8817026 2.69856802 234.5003184 238.0937825
## alpha[2] 247.6981732 2.74071929 242.1202433 245.9699330
## alpha[3] 252.4725536 2.61396199 247.2347767 250.7270517
## alpha[4] 232.6650683 2.58106078 227.3071115 230.9328126
## alpha[5] 231.6676142 2.62784112 226.5553436 229.8802622
## alpha[6] 249.6630372 2.65022329 244.5701236 247.9245898
## alpha[7] 228.7975073 2.83303105 223.1622942 226.7577078
## alpha[8] 248.2318881 2.77037261 242.9606118 246.3422184
## alpha[9] 283.3462482 2.66369313 277.8474352 281.6341046
## alpha[10] 219.1562742 2.68745723 213.8511030 217.4226411
## alpha[11] 258.1861801 2.82267251 252.4833096 256.4499143
## alpha[12] 228.1907674 2.67821382 223.0705725 226.3842512
## alpha[13] 242.3797189 2.79698660 237.1069528 240.4358546
## alpha[14] 268.2588684 2.71732136 262.5688968 266.4373825
## alpha[15] 242.9239134 2.81060606 237.6852887 241.0575444
## alpha[16] 245.3760676 2.72130079 240.1230314 243.5545330
## alpha[17] 232.1473781 2.85795292 226.6159546 230.2257737
## alpha[18] 240.5631123 2.70221334 235.0794508 238.8687811
## alpha[19] 253.7390085 2.68243814 248.3902020 252.0113354
## alpha[20] 241.7204001 2.86416719 236.0308973 239.8677321
## alpha[21] 248.4956916 2.78723439 243.1073384 246.5444261
## alpha[22] 225.2258270 2.81762060 219.5195535 223.3754398
## alpha[23] 228.4431641 2.73461328 222.6051844 226.6854269
## alpha[24] 245.1508235 2.63501708 240.1351554 243.3596567
## alpha[25] 234.5263087 2.78128331 228.8581276 232.7482953
## alpha[26] 253.9806078 2.73452628 248.4030181 252.2410027
## alpha[27] 254.2860958 2.54808777 249.5777908 252.5757730
## alpha[28] 242.9664656 2.64060098 237.5545891 241.3097817
## alpha[29] 217.8840352 2.81530953 212.1043874 216.0266376
## alpha[30] 241.3974418 2.88412281 235.8968942 239.4616713
## beta[1] 6.0701699 0.23988837 5.6295599 5.9027579
## beta[2] 7.0579168 0.25275099 6.5518377 6.8862446
## beta[3] 6.4848320 0.26529857 5.9752814 6.2994094
## beta[4] 5.3533972 0.25641266 4.8134621 5.1852477
## beta[5] 6.5707013 0.24824578 6.0609056 6.4080001
## beta[6] 6.1633262 0.23565976 5.7099808 6.0048494
## beta[7] 5.9698319 0.25569841 5.4560169 5.7957586
## beta[8] 6.4179816 0.25450576 5.9172059 6.2578005
## beta[9] 7.0513725 0.26580336 6.5452452 6.8709382
## beta[10] 5.8596212 0.24138057 5.3636823 5.7065060
## beta[11] 6.8038397 0.25512535 6.2984234 6.6369110
## beta[12] 6.1222221 0.24219197 5.6385177 5.9814465
## beta[13] 6.1698860 0.25444626 5.6576177 6.0118561
## beta[14] 6.6862953 0.25884672 6.1883825 6.5049567
## beta[15] 5.4226856 0.23724832 4.9625396 5.2628071
## beta[16] 5.9273245 0.23856428 5.4530161 5.7643915
## beta[17] 6.2684769 0.24284397 5.7956838 6.1145256
## beta[18] 5.8443937 0.25438384 5.3403493 5.6773952
## beta[19] 6.4129012 0.25713723 5.9076234 6.2469126
## beta[20] 6.0524297 0.23421101 5.6010229 5.8954541
## beta[21] 6.4014545 0.24194834 5.8873367 6.2503140
## beta[22] 5.8727762 0.23951242 5.3987479 5.7130497
## beta[23] 5.7385003 0.25005281 5.2601348 5.5649520
## beta[24] 5.8906368 0.25106056 5.3839170 5.7282164
## beta[25] 6.8999612 0.23893841 6.4272533 6.7421982
## beta[26] 6.5443404 0.24572569 6.0718658 6.3777420
## beta[27] 5.9074945 0.23014974 5.4374604 5.7517532
## beta[28] 5.8535090 0.24414947 5.3956719 5.6839514
## beta[29] 5.6820615 0.24489840 5.2082358 5.5109580
## beta[30] 6.1387519 0.23648680 5.6925340 5.9719287
## mu_alpha 242.5178723 2.68611680 237.2267209 240.7470762
## mu_beta 6.1907515 0.10586261 5.9768566 6.1185673
## sigmasq_y 37.7777138 5.88761902 27.8452809 33.5348247
## sigmasq_alpha 214.5946270 61.49038886 124.9580693 169.8228966
## sigmasq_beta 0.2702446 0.09767191 0.1280761 0.2027839
## sigma_y 6.1281324 0.47321359 5.2768627 5.7909261
## sigma_alpha 14.5076174 2.03169713 11.1784636 13.0316111
## sigma_beta 0.5121260 0.08932866 0.3578772 0.4503154
## alpha0 106.3213395 3.55335578 99.5836386 103.8685244
## lp__ -439.0926770 6.64241703 -452.7594430 -443.3948416
## stats
## parameter 50% 75% 97.5%
## alpha[1] 239.7953230 241.8226917 245.1152715
## alpha[2] 247.7401478 249.5451986 252.9416594
## alpha[3] 252.5223337 254.2400050 257.4981371
## alpha[4] 232.7265422 234.3602142 237.5521243
## alpha[5] 231.5837292 233.3976848 236.8976075
## alpha[6] 249.6935576 251.4262422 255.0387350
## alpha[7] 228.8267203 230.8046293 234.0596074
## alpha[8] 248.1138171 250.0857366 253.9541610
## alpha[9] 283.4004550 285.1599414 288.2149802
## alpha[10] 219.1309693 220.9094533 224.5082897
## alpha[11] 258.1539261 260.1009515 263.6489366
## alpha[12] 228.2154172 229.9946067 233.5722652
## alpha[13] 242.3646350 244.2949197 247.8227405
## alpha[14] 268.3734870 270.1208399 273.3090060
## alpha[15] 242.8941362 244.8184121 248.3579469
## alpha[16] 245.3857135 247.2367111 250.7617351
## alpha[17] 232.2593307 234.0816229 237.5371577
## alpha[18] 240.5258608 242.2339719 246.0926019
## alpha[19] 253.8017039 255.4674728 259.0434919
## alpha[20] 241.8444766 243.6287490 247.2336380
## alpha[21] 248.5045423 250.2830555 254.1742515
## alpha[22] 225.1548806 227.0718694 230.9600325
## alpha[23] 228.4623056 230.2020363 233.9244151
## alpha[24] 245.1735893 246.8749724 250.4368111
## alpha[25] 234.4709175 236.3795426 240.0058128
## alpha[26] 254.0052244 255.7341589 259.6545240
## alpha[27] 254.3213470 255.9263773 259.7263923
## alpha[28] 242.9339396 244.7512727 247.9803880
## alpha[29] 217.8812788 219.6377462 223.5253426
## alpha[30] 241.3461776 243.3199871 247.0525573
## beta[1] 6.0570743 6.2341566 6.5153920
## beta[2] 7.0544968 7.2406474 7.5564443
## beta[3] 6.4859866 6.6679818 6.9935676
## beta[4] 5.3608791 5.5366693 5.8333519
## beta[5] 6.5774503 6.7354097 7.0472251
## beta[6] 6.1621271 6.3159696 6.6412220
## beta[7] 5.9749113 6.1376668 6.4706856
## beta[8] 6.4173170 6.5844322 6.9228151
## beta[9] 7.0489227 7.2211352 7.5893657
## beta[10] 5.8620173 6.0250047 6.3281135
## beta[11] 6.8126315 6.9794234 7.3029219
## beta[12] 6.1214102 6.2774428 6.6014095
## beta[13] 6.1615392 6.3404230 6.6571696
## beta[14] 6.6820962 6.8654986 7.1703971
## beta[15] 5.4274003 5.5803113 5.9101727
## beta[16] 5.9262474 6.0941318 6.3766716
## beta[17] 6.2631793 6.4260049 6.7277244
## beta[18] 5.8406731 6.0096212 6.3422636
## beta[19] 6.4196736 6.5734420 6.9136626
## beta[20] 6.0526775 6.2022182 6.5252676
## beta[21] 6.4127536 6.5604219 6.8741565
## beta[22] 5.8765874 6.0283719 6.3523970
## beta[23] 5.7317950 5.9142991 6.2284916
## beta[24] 5.8945008 6.0516055 6.3634345
## beta[25] 6.8958102 7.0538883 7.3800807
## beta[26] 6.5495454 6.7123946 7.0133975
## beta[27] 5.9087900 6.0611418 6.3470518
## beta[28] 5.8439617 6.0185156 6.3501240
## beta[29] 5.6798621 5.8426826 6.1355987
## beta[30] 6.1407903 6.3061317 6.6028010
## mu_alpha 242.5758159 244.2478887 247.5292893
## mu_beta 6.1941616 6.2561782 6.4073906
## sigmasq_y 37.4916566 41.3133595 51.0177465
## sigmasq_alpha 203.7956678 247.1597949 361.9708566
## sigmasq_beta 0.2524830 0.3159583 0.5210212
## sigma_y 6.1230431 6.4275469 7.1426708
## sigma_alpha 14.2757020 15.7213166 19.0255287
## sigma_beta 0.5024766 0.5621017 0.7218180
## alpha0 106.2670518 108.7476813 113.4326947
## lp__ -439.0995188 -434.6911579 -426.1944556
##
## , , chains = chain:4
##
## stats
## parameter mean sd 2.5% 25%
## alpha[1] 239.8868051 2.64686717 234.5182418 238.2063914
## alpha[2] 247.7819390 2.63128207 242.7569624 245.9904534
## alpha[3] 252.3890607 2.66572282 247.2300563 250.6091994
## alpha[4] 232.5113827 2.74715931 227.2237380 230.5594109
## alpha[5] 231.6458566 2.53605373 226.8579344 229.8967219
## alpha[6] 249.8124634 2.77518481 244.3703152 247.8663916
## alpha[7] 228.7646347 2.54006449 223.8797812 227.0709479
## alpha[8] 248.2667698 2.73245767 243.4232575 246.3657291
## alpha[9] 283.3117646 2.77569442 277.6588566 281.5656057
## alpha[10] 219.2412076 2.59880496 213.9553585 217.4603954
## alpha[11] 258.2921712 2.72258137 253.2055014 256.4278565
## alpha[12] 228.1411716 2.76349840 222.8143739 226.2709706
## alpha[13] 242.3450226 2.72195205 236.9418525 240.5609815
## alpha[14] 268.2169669 2.67339446 263.1405667 266.4143105
## alpha[15] 242.7859829 2.70245770 237.6044974 240.9980768
## alpha[16] 245.2904552 2.60362460 240.0941027 243.5718959
## alpha[17] 232.2875888 2.83258217 226.6480463 230.4321124
## alpha[18] 240.5595852 2.82707994 235.0296508 238.7500433
## alpha[19] 253.7559447 2.54674784 248.7635450 252.0601418
## alpha[20] 241.7166269 2.63303689 236.4060387 240.0098012
## alpha[21] 248.6770200 2.58400847 243.9980593 246.9395452
## alpha[22] 225.2592710 2.56473381 220.0457234 223.5380858
## alpha[23] 228.6543331 2.54869403 223.3204606 227.0181786
## alpha[24] 245.0698896 2.74366813 239.8496764 243.2248450
## alpha[25] 234.3775581 2.69530248 229.0634876 232.6655846
## alpha[26] 253.9478189 2.76844900 248.2238886 252.1426363
## alpha[27] 254.2784326 2.78473953 248.8071095 252.4381398
## alpha[28] 242.9820091 2.83369821 237.1988045 241.0601964
## alpha[29] 217.9590093 2.59007205 213.2464454 216.3117548
## alpha[30] 241.5164540 2.54048044 236.8270761 239.7026616
## beta[1] 6.0580181 0.24627218 5.5857527 5.9008123
## beta[2] 7.0507565 0.25620577 6.5189281 6.8812093
## beta[3] 6.4854125 0.23753143 6.0245502 6.3285150
## beta[4] 5.3398624 0.24777866 4.8290452 5.1826931
## beta[5] 6.5732265 0.23704110 6.0949606 6.4124515
## beta[6] 6.1869150 0.23814487 5.7175528 6.0247618
## beta[7] 5.9707276 0.23424534 5.5030867 5.8181361
## beta[8] 6.4185847 0.24347179 5.9416538 6.2553562
## beta[9] 7.0512614 0.25495670 6.4992986 6.8924669
## beta[10] 5.8553670 0.24150663 5.3953813 5.6984729
## beta[11] 6.7932100 0.26659577 6.2873007 6.6201832
## beta[12] 6.1137323 0.24240055 5.6357968 5.9548262
## beta[13] 6.1788814 0.25917489 5.6451282 5.9986324
## beta[14] 6.6926700 0.25305990 6.2142066 6.5259580
## beta[15] 5.4141822 0.24753669 4.9464824 5.2415233
## beta[16] 5.9195157 0.25222635 5.4302144 5.7419377
## beta[17] 6.2773264 0.26090743 5.7459512 6.1085566
## beta[18] 5.8520853 0.22673299 5.4134053 5.7012770
## beta[19] 6.3978346 0.23774225 5.9293673 6.2328895
## beta[20] 6.0502114 0.25365712 5.5867038 5.8746658
## beta[21] 6.4102521 0.26161467 5.9213165 6.2149224
## beta[22] 5.8473829 0.24024492 5.3953873 5.6748722
## beta[23] 5.7499118 0.26039886 5.2631228 5.5749740
## beta[24] 5.8895225 0.25465043 5.4079135 5.7158520
## beta[25] 6.9072118 0.25400748 6.4352042 6.7312651
## beta[26] 6.5508885 0.24454974 6.0663607 6.3979682
## beta[27] 5.8985094 0.24214781 5.4132212 5.7330473
## beta[28] 5.8462425 0.23848371 5.3683216 5.6908316
## beta[29] 5.6739886 0.24597519 5.1870620 5.5182602
## beta[30] 6.1333733 0.23854329 5.6825510 5.9628392
## mu_alpha 242.4973392 2.66725037 237.2370232 240.7619799
## mu_beta 6.1832235 0.11038014 5.9621103 6.1059631
## sigmasq_y 37.2777193 5.95125437 26.8089044 33.2285181
## sigmasq_alpha 214.8608747 60.80616835 124.6977843 171.5266041
## sigmasq_beta 0.2705845 0.10141391 0.1238200 0.2000185
## sigma_y 6.0867860 0.47852275 5.1777316 5.7644183
## sigma_alpha 14.5208361 2.00254961 11.1668151 13.0968152
## sigma_beta 0.5117773 0.09315135 0.3518807 0.4472342
## alpha0 106.4664232 3.64097096 99.1693945 103.9998035
## lp__ -438.4282561 7.41528386 -453.3311272 -442.8941137
## stats
## parameter 50% 75% 97.5%
## alpha[1] 239.9174372 241.5845400 244.8524258
## alpha[2] 247.6560853 249.6002015 253.1645567
## alpha[3] 252.3512595 254.1966649 257.8076243
## alpha[4] 232.4820886 234.4275281 237.7855349
## alpha[5] 231.5443774 233.3268598 236.5073489
## alpha[6] 249.8912492 251.8065450 255.0849111
## alpha[7] 228.7433729 230.3790644 233.8858232
## alpha[8] 248.2352212 250.1419205 253.7210835
## alpha[9] 283.4241544 285.1674916 288.4733256
## alpha[10] 219.2943985 220.9975449 224.2669334
## alpha[11] 258.2953720 260.0198678 263.8958646
## alpha[12] 228.2183191 230.0537923 233.5704259
## alpha[13] 242.3452957 244.1968484 247.5040710
## alpha[14] 268.2800079 269.9695423 273.5082773
## alpha[15] 242.7711752 244.5444550 247.9377072
## alpha[16] 245.2804796 246.9843015 250.3325518
## alpha[17] 232.2811229 234.0948592 237.8684299
## alpha[18] 240.5480460 242.3804493 246.2023281
## alpha[19] 253.7197835 255.4343062 258.8152436
## alpha[20] 241.7286483 243.4321528 246.7811340
## alpha[21] 248.6886099 250.4149197 253.8052086
## alpha[22] 225.3348428 226.8534969 230.2081357
## alpha[23] 228.7477778 230.3186177 233.5241400
## alpha[24] 245.1033207 247.0040683 250.3501972
## alpha[25] 234.3565560 236.0984075 239.5937414
## alpha[26] 253.9691438 255.7616071 259.3976438
## alpha[27] 254.2556889 256.1721080 259.5763990
## alpha[28] 242.8898041 244.8974316 248.4319710
## alpha[29] 217.8123672 219.6836418 223.1868211
## alpha[30] 241.5028561 243.2711098 246.2570559
## beta[1] 6.0515695 6.2163460 6.5399459
## beta[2] 7.0570363 7.2293432 7.5131586
## beta[3] 6.4873006 6.6485816 6.9524299
## beta[4] 5.3342232 5.4935655 5.8489087
## beta[5] 6.5720259 6.7352865 7.0303646
## beta[6] 6.1807716 6.3455350 6.6586552
## beta[7] 5.9737921 6.1307786 6.4177977
## beta[8] 6.4112710 6.5856315 6.8964935
## beta[9] 7.0448015 7.2179397 7.5485346
## beta[10] 5.8528694 6.0124665 6.3158228
## beta[11] 6.8005100 6.9565755 7.3361636
## beta[12] 6.1107271 6.2756795 6.6109746
## beta[13] 6.1763173 6.3518566 6.6787559
## beta[14] 6.6836511 6.8541864 7.2174820
## beta[15] 5.4187277 5.5824149 5.9140072
## beta[16] 5.9254710 6.1021868 6.3734574
## beta[17] 6.2892993 6.4540046 6.7626293
## beta[18] 5.8488151 5.9985340 6.3237274
## beta[19] 6.4052169 6.5687162 6.8381995
## beta[20] 6.0458495 6.2312275 6.5549838
## beta[21] 6.4128112 6.5982596 6.8794340
## beta[22] 5.8572477 6.0135790 6.3153987
## beta[23] 5.7372930 5.9283252 6.2819208
## beta[24] 5.8972562 6.0478470 6.4172307
## beta[25] 6.9033780 7.0790977 7.4214783
## beta[26] 6.5489587 6.7149183 7.0411321
## beta[27] 5.8986770 6.0590043 6.3676260
## beta[28] 5.8435118 6.0138741 6.3251701
## beta[29] 5.6840831 5.8287323 6.1522461
## beta[30] 6.1343568 6.2967186 6.5967581
## mu_alpha 242.5189834 244.3516290 247.6480769
## mu_beta 6.1864496 6.2580496 6.3874073
## sigmasq_y 36.7048170 41.0566691 49.5044001
## sigmasq_alpha 206.6644391 247.8234694 362.3880683
## sigmasq_beta 0.2517077 0.3223430 0.5208619
## sigma_y 6.0584500 6.4075473 7.0359363
## sigma_alpha 14.3758282 15.7424099 19.0364910
## sigma_beta 0.5017048 0.5677526 0.7217076
## alpha0 106.5240975 108.9072703 113.9302470
## lp__ -438.0246915 -433.3851105 -425.2839243
# stan script
eightschools<-"
// saved as 8schools.stan
data {
int<lower=0> J; // number of schools
real y[J]; // estimated treatment effects
real<lower=0> sigma[J]; // standard error of effect estimates
}
parameters {
real mu; // population treatment effect
real<lower=0> tau; // standard deviation in treatment effects
vector[J] eta; // unscaled deviation from mu by school
}
transformed parameters {
vector[J] theta = mu + tau * eta; // school treatment effects
}
model {
target += normal_lpdf(eta | 0, 1); // prior log-density
target += normal_lpdf(y | theta, sigma); // log-likelihood
}"
schools_dat <- list(J = 8,
y = c(28, 8, -3, 7, -1, 1, 18, 12),
sigma = c(15, 10, 16, 11, 9, 11, 10, 18))
fit <- rstan::stan(model_code = eightschools,
data = schools_dat)
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
print(fit)
## Inference for Stan model: 432a5acc93a981154b518b22530ee9a6.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 7.84 0.11 5.01 -2.32 4.58 7.87 11.14 17.95 2145 1.00
## tau 6.63 0.14 5.52 0.34 2.50 5.26 9.24 21.06 1568 1.00
## eta[1] 0.42 0.02 0.95 -1.54 -0.21 0.45 1.04 2.21 3742 1.00
## eta[2] -0.03 0.01 0.86 -1.74 -0.59 -0.03 0.52 1.65 3510 1.00
## eta[3] -0.18 0.02 0.91 -1.98 -0.79 -0.20 0.42 1.64 3361 1.00
## eta[4] -0.02 0.01 0.88 -1.76 -0.58 -0.03 0.53 1.73 3935 1.00
## eta[5] -0.35 0.02 0.87 -2.04 -0.92 -0.38 0.23 1.37 3198 1.00
## eta[6] -0.22 0.01 0.89 -1.93 -0.81 -0.22 0.36 1.56 3998 1.00
## eta[7] 0.35 0.01 0.85 -1.37 -0.18 0.37 0.91 2.04 3681 1.00
## eta[8] 0.08 0.02 0.94 -1.77 -0.55 0.09 0.71 1.93 3914 1.00
## theta[1] 11.57 0.16 8.44 -2.06 6.21 10.32 15.83 31.71 2920 1.00
## theta[2] 7.79 0.09 6.26 -4.37 3.88 7.65 11.68 20.55 4550 1.00
## theta[3] 6.16 0.13 7.67 -10.65 2.04 6.59 10.70 21.09 3349 1.00
## theta[4] 7.60 0.10 6.44 -5.40 3.63 7.60 11.62 20.48 3772 1.00
## theta[5] 5.05 0.11 6.27 -8.44 1.21 5.50 9.28 16.06 3451 1.00
## theta[6] 6.16 0.10 6.65 -8.22 2.28 6.61 10.40 18.41 4175 1.00
## theta[7] 10.68 0.12 6.68 -0.78 6.22 10.03 14.55 25.71 3366 1.00
## theta[8] 8.51 0.14 7.91 -7.77 3.98 8.27 12.68 25.18 3414 1.00
## lp__ -39.41 0.07 2.61 -45.24 -40.96 -39.17 -37.59 -34.99 1267 1.01
##
## Samples were drawn using NUTS(diag_e) at Mon Jan 7 16:18:43 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
plot(fit)
## 'pars' not specified. Showing first 10 parameters by default.
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)
pairs(fit, pars = c("mu", "tau", "lp__"))
la <- extract(fit, permuted = TRUE) # return a list of arrays
mu <- la$mu
### return an array of three dimensions: iterations, chains, parameters
a <- extract(fit, permuted = FALSE)
### use S3 functions on stanfit objects
a2 <- as.array(fit)
m <- as.matrix(fit)
d <- as.data.frame(fit)
http://www.stephenpettigrew.com/r/
https://cran.r-project.org/web/packages/apaTables/vignettes/apaTables.html
papaja: Reproducible APA manuscripts with R Markdown